{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": "true"
   },
   "source": [
    "# Saddle plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:50:23.743811Z",
     "start_time": "2018-09-25T16:50:22.357663Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "import bioframe\n",
    "\n",
    "import cooler\n",
    "\n",
    "import cooltools\n",
    "import cooltools.eigdecomp\n",
    "import cooltools.expected\n",
    "import cooltools.saddle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:51:12.838618Z",
     "start_time": "2018-09-25T16:51:09.341801Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "--2019-10-30 15:16:16--  ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE93nnn/GSE93431/suppl/GSE93431_NIPBL.200kb.cool.HDF5.gz\n",
      "           => ‘/tmp/GSE93431_NIPBL.200kb.cool.gz’\n",
      "Resolving ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)... 130.14.250.7, 2607:f220:41e:250::13\n",
      "Connecting to ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)|130.14.250.7|:21... connected.\n",
      "Logging in as anonymous ... Logged in!\n",
      "==> SYST ... done.    ==> PWD ... done.\n",
      "==> TYPE I ... done.  ==> CWD (1) /geo/series/GSE93nnn/GSE93431/suppl ... done.\n",
      "==> SIZE GSE93431_NIPBL.200kb.cool.HDF5.gz ... 37884973\n",
      "==> PASV ... done.    ==> RETR GSE93431_NIPBL.200kb.cool.HDF5.gz ... done.\n",
      "Length: 37884973 (36M) (unauthoritative)\n",
      "\n",
      "GSE93431_NIPBL.200k 100%[===================>]  36.13M  14.8MB/s    in 2.4s    \n",
      "\n",
      "2019-10-30 15:16:19 (14.8 MB/s) - ‘/tmp/GSE93431_NIPBL.200kb.cool.gz’ saved [37884973]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# download a Hi-C dataset from Schwarzer et.al. \"Two independent modes of chromosome organization are revealed by cohesin removal\", 2017\n",
    "\n",
    "!wget ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE93nnn/GSE93431/suppl/GSE93431_NIPBL.200kb.cool.HDF5.gz -O /tmp/GSE93431_NIPBL.200kb.cool.gz\n",
    "!gunzip -f /tmp/GSE93431_NIPBL.200kb.cool.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:51:13.865263Z",
     "start_time": "2018-09-25T16:51:13.842147Z"
    }
   },
   "outputs": [],
   "source": [
    "coolpath = '/tmp/GSE93431_NIPBL.200kb.cool'\n",
    "c = cooler.Cooler(coolpath)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:51:52.903044Z",
     "start_time": "2018-09-25T16:51:52.891798Z"
    }
   },
   "outputs": [],
   "source": [
    "# Define continuous genomic regions for calculations of contact frequency decay \n",
    "# with distance (aka \"expected\").\n",
    "# Typically, we calculate expected separately for each chromosomal arm because\n",
    "# centromeres additionally suppress contacts accross themselves.\n",
    "# In mice, chromosomes are acrocentric and expected can be calculated \n",
    "# for whole chromosomes.\n",
    "\n",
    "regions = [(chrom, 0, c.chromsizes[chrom]) for chrom in c.chromnames]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:52:19.033257Z",
     "start_time": "2018-09-25T16:52:01.255106Z"
    }
   },
   "outputs": [],
   "source": [
    "# Download and compute gene count per genomic bin \n",
    "\n",
    "bins = c.bins()[:]\n",
    "genecov = bioframe.tools.frac_gene_coverage(bins, 'mm9')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:57:03.355125Z",
     "start_time": "2018-09-25T16:52:20.976856Z"
    }
   },
   "outputs": [],
   "source": [
    "# Perform eigenvector decomposition in cis, sorting and flipping eigenvectors \n",
    "# according to their correlation with the number of genes in each bin.\n",
    "\n",
    "cis_eigs = cooltools.eigdecomp.cooler_cis_eig(\n",
    "    c, \n",
    "    genecov, \n",
    "    regions=None, \n",
    "    n_eigs=5, \n",
    "    phasing_track_col='gene_count')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:57:03.511710Z",
     "start_time": "2018-09-25T16:57:03.357213Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'chr1 position, bp')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3cAAACaCAYAAAAKAbHXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9d3xb9b3//zra1rI1LMl7xHYSx9mEhCxmKGEVSmkpFDrg293S8bulvfd23fbSS28p0Ht7u5gttBQKLRASRoAkZE/HGR7xkLxla1p7nt8f0lE8JFuSj5bzeT4eeSTR/NjSOefzXq8XRdM0CAQCgUAgEAgEAoFQ2HByvQACgUAgEAgEAoFAIMwfEtwRCAQCgUAgEAgEwgKABHcEAoFAIBAIBAKBsAAgwR2BQCAQCAQCgUAgLABIcEcgEAgEAoFAIBAICwAS3BEIBAKBQCAQCATCAoCX6wWkglqtpmtra3O9DAKBQCAQCAQCgUDICSdOnDDRNF0a776CCu5qa2tx/PjxXC+DQCAQCAQCgUAgEHICRVGGRPeRtkwCgUAgEAgEAoFAWACQ4I5AIBAIBAKBQCAQFgAkuCMQCAQCIQd0jE7gmy+ewqsnB+H0BXO9HAKBQCAsAApq5o5AIBAIhIXCm20j+GfrMP7ZOgwh7wyuX6bDj25phloqzPXSCAQCgVCgkModgUAgEAg5QG92o0pZhFe+fAXuWleFd86N4rbfHEDnqCPXSyMQCARCgUKCOwKBQCAQckC/2YValQRra5T4yUdb8NIXr4A/GMYdvz2IXWdG0DnqwAmDFaf6raBpOtfLJRAIBEIBQNoyCQQCgUDIAXqzG7esLIv9f2VVCV772iY88NxxfPmFk1Me+4+vbMTqakW2l0ggEAiEAoMEdwQCgUAgZBmb2w+7J4AapWTK7WXFRXj5S1dgd/sYuBQFbyCE77x8Gl1GBwnuCAQCgTAnJLgjEAgEAiHLGMxuAECNSjzjPrGAh1tXlgMAgqEwvvdqG/pM7qyuj0AgEAiFCZm5IxAIBAIhyxgsTHAnmfVxPC4HVUox+kzObCwrZZ4/bMBHf3OAzAQSCARCnkCCOwKBQCAQsozB5AIAVCtnVu6mU6+WQJ+nlbtdZ0dwesAG44Qv10shEAiXMMYJL+yeAEk0gbRlzhu7JwBfMASNTJTrpRAIBAKhQDBY3NDKhSgScOd8bK1Kgg8vmBAO0+BwqCysLjlCYRqt/TYAQPvoBHTF5DpIIBCyT9ugDbf+7wEAgFjARbVSjIc/thxrLtE5ZVK5myeff/YYvvW31lwvg0AgEAgFhMHsmrMlk6FWLYEvGMbohDfDq0qNLqMDLn8IANAxQrz5CARCbvjwggkA8N0bFuOuddVweIP4wp+OY9Canx0PmYYEd/NEIRbA6grkehkEAoFAKCAMZjdqkmjJBCJtmQDQF23lzBdO9lsBAEV8LjpGJ3K8GgKBcKlytM+CJq0UX7mqAT+8pRnPff5y+IJhPPDccbh8wVwvL+uQ4G6eKCV8WN3+XC+DQCAQCAWC2x/EmMMXVykzHrX5GtwZbFBJBLhikQrtIyS4IxAI2ScUpnHCYMW6WmXstgaNFL+5ew26jA48+GIrrC7/JTWLR2bu5olCLIAl+qWhqPyZhSAQCARCftKfpFImg04ugpDHgX6ewd1/vnkeGxepcfUSzbxeh+FUvxWrqxVYrJNib9c4fMEQhLy5ZwgJBAKBLdpHJuD0BXF5nXLK7VubSvHDm5vx4zfOY/VP34WAx4FWLoREwIOQz4WIx8F1S7X49IaapGafC4mcBncURT0N4GYAYzRNt+RyLemikAjgC4bhCYQgFpBYmUAgEAizwyhfJlu543Ao1Kkl86rc2T0B/PHDPvSOu1gJ7qwuP3pNLnz8skpUKcQIhWl0jzmxrLx43q9NIBAIyXK0zwIAUyp3DJ/ZWIsGjQydRgeME14YJ7zw+EPwBcOwuPz4z53t+P2+Xnz5qkX45LoqSIUX9/EXjA48faAPdWoJvrB1UdZ+HjbIdTTyLID/BfCnHK8jbZRiAQDA6g4UfHDnC4YwNuGDRi4k2VcCgUDIEP2WSJBWo0yucgdEFDO7xtIXLTk3bAcAnB60sdJpcmogMm+3ploBtVQIICKqQoI7AoGQTY7pLagoKUJ5SdGM+yiKwuZGNTY3qhM+97F3u/DTHefx853tWFOjwKZFapzst2Jv1ziEPA7+35b6TP8IrJPTaISm6X0URdXmcg3zRSGJBncuPyrifLEKgR/88yzePDMCiysyO3j76go89slVOV4VgUAgLEz0ZjdKxHwUi/lJP6dWLcF7HUYEQ2HwuKmPy58biszEmZx+DFo9qEpSzCURJw02cDkUVlQWQ8DlQMjjkLk7AqFA8PhD+NW7nfjKVQ2xfWwhQtM0jukt2NpYmtbz19Uq8Zf/twEnDFbsbjdib+c4HtvdBbVUiO9sa8I9G2qgLMDfT2GXmvIARfTizARGhchrrUOoUorx2Y21eLNtBN1jzlwviUAgEBYs/WZ30vN2DHVqMQIhGkM2T8rPBYAzQ3ZwOVTEm27ANv/grt+KpWWyWMdKk1aGjlFih0AgFAI7z4zgjx/2YVl5MW5bXZHr5aRNn8kFk9OPdXUzWzJTYW2NAmtrFHjohiWwuf0QC3gQ8ApXczLvV05R1BcoijpOUdTx8fHxXC9nBrHKXYEqZrp8QUx4g7hpRRm+cW0jVlWV5J2XEoFAICwk9GZX0jYIDHVqKYD0FTPPDtuxtVENIY+D1gFbWq/BEArTOD1gm2IQvEQnI3YIBEKBsOvsCABgzFHY+71j+sTzdulSIhYUdGAHFEBwR9P0H2iavoym6ctKS9Mru2YSZuauUCt3TCBXViwCAGiLRTA5fQiEwrlcFoFAIGQNuzuAYJbOef5gGMM2D2qTFFNhqFVHHp+OYqbTF0SfyYVVVQq0VBTj9DyDu87RiHn55OBuaZkcJqcf4w7fvF6bQCBkFoc3gH1dEdPvsYnCPl6P9lmhkgiwqDT1boaFTN4Hd/mOvIgPDhURVClERmxMcFcU/VsEmgbGyAWaQCBcAgRCYVz1yw/w7EF9Vt5v0OpGmAaqU2ytLJUKIRFw06rcnR+eAE0DLRVyrKoqwZkh+7wSeIx5+erqkthtS8pkAECqdwRCnvNe+xj8oTC4HKrg93rH9BZcVqsgVmTTyGlwR1HUXwEcArCYoqhBiqLuz+V60oHLoVBcxIe1QCt3I3YPgIuVO5088veovbBL9QQCgZAMQ1YPrO4Azg1nJygxRD3uUq3cURSFulIJ+szulN/z7FBEKXN5RTFWVpXAFwyjcx7zcXs6x1AqE6J6UmvpEp0cAIioCuGS4ILRUbAdTjvPjEAnF2FVVUlBV9qNE170W9ystmQuFHIa3NE0/SmapstomubTNF1J0/RTuVxPuigkAlgKdOaOCeK00aBOS4I7AoFwCdFnjlTCGGPxTGOIVt6qUwzugIgdQjptmWeH7CiVCaGRi7C6KlJtO5Vma2afyYX3OsZw17qqKdlypUQArVyIjhEiqkJY2JicPmx/4kP87dhArpeSMk5fEHu6xnFDiw5aubBgZ+5omsbju7sAAJsa4tscXMqQtkwWUIoFhVu5m/BCJRFAxI/42jEVPCKqQiAQLgWYYMmQRkUsHXrGXZAKeSiNesOlQp1agkGrG/5gahWDs8N2tJRHKmuViiKoJIK05+6e3t8HPoeDe6+omXHf0jI52oliJmGB0zvuQjBMx7wjC4n3O8bgD4Zx04oyaGSigm3LfP6wAX89OoCvXr0IS8vkuV5O3kGCOxZQSAQFPHPngS4a0AFAiZgPAY8DIwnuCATCJQAT1JmcPrj9wYy/X5fRgSatNK0ZkTq1BGE6tSqjxx9C95gTyysi5uIURWFlVUlaiplWlx8vnxjAR1eVQyMTzbh/aZkc3WMOGMxTq4u940489m5XykEpgZCPMN/vC8bCs43a2TYCjUyItdUKlMqEcHiD8AZCuV5WShzqMeMnb5zHtUs0+M62xbleTl5CgjsWUIgLeebOG6vWAZELv04uwghpyyQQCJcAkwVKMt2aSdM0uowOLNbJ0np+oybyPEbQJBnOj0wgTAPLosEdAKyqKkHPuBMT3tSSkn852g9vIIwHttTHvf/TG2pQxOfi6389FQvkrC4/PvvMMTzx3gXsaBtO6f0IhHyESQhdGHOCpukcryZ5XL4gPugcw/YWHTgcCqWySPdAISlmjtq9+MoLJ1CrluDxu1aBwyFCKvEgwR0LMDN3hXSQM4xOeKdU7gBAVyyCkQR3CwKLy0/mJwmEWdCbXTFhkEy3Zo47fbC6A7EgLVVaKuRYVCrBC4cNST+HaR1rmRbc0TTQNpB8W5kvGMKzB/XY2lSaMDitKCnCLz6+Am2DdvzirQ4EQmF85YWTGLV7oZOL8Nyh5NdNIOQrjCiS3RPAuLNwAqMPOsfgC4axfXkZAEDDBHcFNHf3XocRVncAv75rNWQifq6Xk7eQ4I4FlGIB/MEw3P7CKm17/CHY3IGYDQKDTi4iM3cLhH//5xnc8+ThXC+DQMhLAqEwBq0ebG2KDOQPZLhyx7RxpVu5oygKn9lYi9OD9qTbKs8O2aGUCFA+KYm3sjIiqtI6kHwF8I3TIxh3+PDA5rpZH3dDSxnuu6IGT+7vw31PHcWhXjN+/rHl+PJVi3B6wDZvA/VCYdzhwwPPHS9YD1xCYgxmF4qiOgXdBdSauevMKNRSYUxdkmmtLqS5u0GrBzwOlfY59FKBBHcsoJAUppH5dBsEBl1xJLgrxEokYSodIw70jLvS8sYiEBY6g1YPQmEaKytLIBPxMl65Y+wHmrTpb0w+tqYSUiEPf5rky2cwu7DtV3txsMc04/FnhiawrFw+ZcavWMxHrUqMs0PJ2xa8c24U1UoxtjTOrUz3rzcuRXOZHId6zfji1nrcsbYSd6yNrPu5LPkJ5pq9XePY3W6ct2E8If8wmN2xhNCFscII7jz+EN7vGMMNLVpwo62MGnmkcldIdgiDVg/KS4piPwMhPiS4YwGFOBLc2QpMVIVp15vRlikXwR8MF6xIDCFCIBSOzRDt7RzL8WoIhPyDUcqsU0tQoxJnfOauy+iAUiKAWipI+zWkQh7uWFOBHW0jMDl98PhD+OKfT+DCmBPvt089zr2BEC4YHVNaMhmatDJcGEte2XLQ6kGDJjkhGBGfiz9+5jL89KPL8N0blsTW/fG1ldjRNlxQm8l06Yj6/RVS2x5hbmxuP+yeAC6rUaK4iI8uY2Gow+7pHIMnEMKN0ZZMINJ1xuNQBdWWOWR1o1JRNPcDL3FIcMcCSkmk77fQvO4Y0ZTy6W2ZxcTrbiEwaPUgGI5UX/d0jed4NQRC/qGPqt7VqiWoUUoyHtx1Gh1oTDJAmo17r6iFPxTGi0f78a//OINOowNqqQBtQ1Nn6M4NTyAYrUxOp1Erhd7shi+Y3DjBkM2DipLkN1UVJUW494raKRn2+66oQSBE469H+5N+nUKlM7rpNzsLa19AmB2mul+jEqNRIy2Yyt2bZ0agkghw+STDbw6HgloqLChBlUGrhwR3SUCCOxZgKneFppjJtGVOr9zFjMwnPFlfE4E9+kyRi86KymIc7jUXnNwxgZBp9CYXZEIeVBIBqpRiDFrdCIUz045O0zQuGJ2szIo0aKTY3KDGr9/vxj9ODeFb1zXhpuVlODdkn7J+piVwdfXM4K5JK0MoTENvmjugdfqCsHsCqJjnpqq+VIqtTaV44YgBgdDCtkVoH2GCu8LZOBPmZnJCqFErRXcBBHfeQKQl8/plOvC4U7f9pTJhwczceQMhjDl8qFSIc72UvIcEdyygLNiZOy8UYn7MwJwhZmRuL4wDnhCf3vHIRehzm2rhDYRxpM+S4xURCPlFn9mNGrUYFEWhRiVGIETHkl5sM2z3wukLzmvebjL3XVEDfzCMa5do8LWrG7C8sgQufyiW1AGA04M26OSiWMJuMoxiZzJtZUPWyO8klcpdIu5cWwnjhA/tI8nP+xUaJqcPpmhQZy6wfQFhdvqjlbtqpRgNGhksLn/eB/B7Osfh9odw06SWTAZNAQV3w7bIeYhU7uaGBHcsIBfxwaEivdiFxKjdC13xzIOkVCYERYEoZhY4fSYXSsR8bG8pg5DHwd5O0ppJIEzGYHahViUBgJgdQqZaM7tYEFOZzLZmLf5432V44lOrweFQMZPytsGLrZmtAzasrJo5bwcA9aUScKjkBCGGbJHfyXwrd0CknQ0AjAXUCpYqjHAOh0IsyCMsDPRmN3RyEUR8Lho1UgD5L6qy6+wIFGI+NtQrZ9ynkQsxXiAzd4MsJpkWOiS4YwEOh0KJWFCQM3flxTMzunwuB6VSIUYzlMEmZIc+kwt1aglEfC7W16uwp4uIqhAIDIwNwozgLkOKmUyFrEkrZeX1KIrCtmYtpEIeAGBRqQRFfG4suLO6/DCY3VhVpYj7fBGfixqVBBeyXLljqojGBZw8ZKqSKypLyMzdAsNgdqE6mqBo1KYX3H3npdN49J1O1tcWD28ghPfax/CROC2ZAFAqE8Hs8iNYAG3STHBXqSRtmXNBgjuWUIj5sLoKS11yxO6ZMW/HELFDIBnHQoYJ7gDgqqZS9I67Mu7jRSAUCgOWyHxdbfQYKSsWgcehYgbFbNNpdEArF6JEnL5S5mzwuBwsK5fjbFRU5fRgZN4uUeUOiMzuJbMxHbR5IIgm/eaLSiIARRWWt1aqdI46oJYKsVgrI5W7BYbB4kZtNLjTyUWQCnnoTkEx0zjhxaunBrHr7GimljiFt86OwukLxozLp6ORCUHThdE+PGh1g8ehoJXN/zy00CHBHUsoJYKCmrnzBkKwugMzPO4YtHIRjEQts2Bx+4MYsXtRzwR3i0sBENVMAoGBUb2rU0c2ajwuB5WKosy1ZRodrLVkJmJ5ZXFEITMUxukBOygKsXbNeDRppdCbXPAHZ8/aD1k9KCsRgcOCtxSPy4kq9C3c60vHqANLy2RQSSP7gnCGRHoI2cXlC2Lc4UNNtNpPURQaNFJ0pWBk/mbbCGg6knxNVqk2XXzBEB59txNLy+TY0hDfn1ITDZQKQTFzMHoeileBJEyF/IZYokQsgLWA2jIZm4OyODN3kdtFBTNz5wuGcPcfD+Pl4wO5XkrewCjg1aml0b8lqFIWEb87AiFKX9TjjtmoAUCVUpyRtsxQmEb3mDPjwd2KymJ4AiH0jLvQOmBFo0YKmYif8PFNWhmCYTr2u0hEqjYIc6GVCxdsW2YoTKPL6MBirQwqqRDBMI0Jb2F19RDiwyR+mLlRIJIgSaUtc0fbMDjUxXNCJnn+cD8GLB58f/uShIkZTbRNuhC87oZsHlSWkJbMZCDBHUsoCyy4G4kFd4krd3ZPAB5/fsnnf3hhHI5pF8pnD+hxsMeMx97tWvDy2snSN8mcGYhkGNfXqaaILRAIlzJ680UbBIZMGZkPWNzwBsJYnOnKXUXE8uD0oA2nB+1x/e0m0xAThJi9rWyY7eBOJlqwgip6swu+YBhLyuQxs3oTmbtbEBgYG4RJCaFGTaT1NhkrrAGLGyf7bbhtVQWAi8I7mcDuCeB/3r+ALY1qbG0qTfi4UqZyVwBt0oPEwDxpSHDHEgqJAFZXADRdGO0XiTzuGHQxr7v8yeb0m92496mj+PLzJ2NeTianD//7fjcqSoowbPdi55mRHK8yP2Dk0GvVF7NcGpkQFpe/YL6jBEIm0ZvdqFVLphiKVyvFsHsCsLvZrbQwhtaNLImpJKJeLYFEwMWuMyOwuPxYWTV7cLeoVAoOhVnbyvzBMMYcPlaUMhk08sKRX0+Vjqi/3RKdDOrojGK+S+UTkoNp5a6eVLlriB7T3eNzV+HejO5Pvn5tIwRcTkaDu9/t7YHdE8BDNyyZ9XHMHG2+t2X6giEYJ4jHXbKQ4I4llBI+/KEwXHlW6UrESBJtmcDF9s184HCfGQCwv9uEX0aVph57twueQAjPfX4d6tUSPLW/jwQvAHpNLpQViyAW8GK3KSWCSIuQJ5jDlREWGv5guCCPOb3JNaW9CgCqlZGMPNvVuwux4C6zlTsOh8KyiuLYbO2qOYI7EZ+LaqUY3bNU7kbsHtA0u/LjGpkIZpdvQXZadI5OgENFqqIqUrlbUOjNbiglAsgntTozdgiHe8xzPv+N08NYVVWCOrUEizRSdGQouBuxe/D0/j7ctqoCLbPM3AKAgMeBQszP+7bMYVtkfaRylxwkuGMJRgEtmdJ8PjBq96JEzEeRgBv3fi0T3E3kjx3C0T4LFGI+PnV5FX67pwf/+/4F/PVoPz69oQYNGhk+v7kObYN2HCVm3VOUMhmYLLLJld8ZOkLh4PAGsO2xvbFkS6Hg9gcxaHXPOEYYOwSDZfYZtFTpNDpRqSiK2RZkkhUVxaBpQMjjYLFu7mCyUSubtXIXs0FgcVOllYtA0wvTA6591IH6UilEfC5UkmjljpxzFwT9lpkJoYqSImxcpMKj73bh57vaY11FQMRuhUl89Y47cW54AresLAcALNXJMla5e79jDL5gGF+9uiGpx2tkoryvpA9aIwk3EtwlR+avNJcISia4c/tRVQAeHCN2T6z1Mh6xtkx7/hzwR/ssWFerxI9vXYbzIw788p0uFBfx8eC1jQCAO9ZU4tF3OvHk/j6sr1fleLW5pc/kwk3TpI+V0dkii8uPRYlb8AmEpPmvXR0wmN0pqcXlAy8c7keYvqgiy8C0WxlYFlXpGnVkfN6OYXllJFPfUlEMfhKqck1aKT7oGIM/GIaAN/Pxg7aotxSLQgZaeSToMU74EnaPFCodoxNYEZ11VIj5oChSuVso6E1urKud6htJURSe/dzl+I8d5/D7vb04PzyB1VUlONxrQeuADSqpANtbyuD0BUBRiF2XF+tkePXUEGxuP+v2KAazGwIeJ6aWPRcRI/P82evFg3jcpQap3LGEYtLGuRAYsXsTiqkAgETIg0zEyxtFs1G7F/0WNy6vU0LI4+J3n16DJq0U/3bT0tjvvkjAxac31GB3u3FO9beFjNXlh80dmFGVYFqEyPwHgQ0O95rxwpF+AIX1nXL7g/jd3h5saVRjbY1yyn1SIQ9qqYBVxcxAKIxekzPjLZkMTGAxV0smQ6MmopipN8c/Zw5ZPaCoxPPZ6aCRRRX68uT6whZOXxADFg+WRiumPC4HSrGgoI4PQnx8wRBG7J4p6roMAh4HP7ttOR65YzmO9Frwvx90wxcM4d4rarCsXI7nDxvw0vFBrKtVxo4jpqqeidZMg9mFKkVR0tYlpbJCCO6Ix10qkModSzBVkUJRzByxe2ObgETo5KKY8EquOaqPtFqur4tU5MqKi/DOt66c8bh7r6jB7/b24MVj/fj+9qVZXWO+0BsNbOtL47dlFoJZKSG/8fhD+N4rbahWirFYJ0NXCia+uebPhwwwu/z45nWNce+vUUkSBjrpoDe5EAjRWKzLrJgKQ61KjG9va8JNK+KbFk+HEXlJ5MM3ZPNAKxPFreqlS6xyl+cbylRh2uyW6OSx21RSAcxZqNwxytaJRi0I86NjxIEwPfO6OplPrqvGtmYdeFxqylzehDeAfV3jaC67+L1gviOdow5sYLnTyGB2xw1CE8EEdzRNTxGYyieIx11qkN8SSyjEkQPZ4sp/P5sJbwAWl39G7/h0dMUijOaJgtLRPjOkQh6Wls2e/dbIRKhVSWAwZcaIuBC4aIMwdTOpEDOVu/wP7hzeQMYNXgnp8/juLujNbvzXHctRqSgqiO8UEKna/X5fb9yqHUONSsxqWyajlJlpjzsGiqLwjWsbsag0uWCSUcy8kKC1dsjqQXkJe1U7AFBJheBQC69y1xP1LZusiqqSCLMyc/fgi6fwrb+1Zvx9LlXeOD0MPpfCVU2aWR83XXAFAOQiPm5eUY76ScekVi5EcRGf9codTdPot7hj88PJoJGJ4A+FYWNZJZhNhqzs2rEsdEhwxxJyER8cCrAVQOWOaTmqnSu4k4tgzBO1zKN9FqypUSSVtdHKRTDmufJTJukzOcHjUDMGjwU8DuQiXtZbh72BEG77zQHsv2Ca87FDNg9+8M+zWPuz3fjpjvNZWB0hHZ4/bMCtK8uxcZEaaqkQTl8Q3kD+B+N/OmSAxeXHt7Y1JXxMrUqC0Qkvax6fXaMOcCgkHWxlGxGfixqVBOeG43tgDtk8qGBZfpzLoaCWLjwjc73ZBR6HmrIJVUkFWZm56x53omN0IuPvcykSDtPY0TaCK5tKUSzmz/2EJKAoCot1MnRO+sx++NpZfPul+QXoJqcfbn9ozuT9ZDQF4HU3aPUQG4QUSDu4oyjqD2wupNDhcCgoxIKCmLljWo7mKtvrikUYc3gRzLFctcXlR5fRifV18TPt09HIhHnv2ZJJ+kwuVCvFccUUVFJh1hXqzgzZ0Tpgw+5246yPe2L3BVz5iw/w4rF+FPG5OD9MNir5iNsfhMsfQnN5pK3oolFzfh9zHn8If9jXiyubSrGmWpHwcbVqdu0QuoxO1KokEPHzt13uyqZS7LtggtM31SYlHKYxYs9Mxlwrz3+FvlTRmyPn3slJSHWWzrkmhw/Ddi/C4cKzJcl3juotGJ3wxpQu2WKJLqJUS9M0This+NMhA949b5yXtUy/hdnfpRPc5WeyxRcMwejwEqXMFJg1uKMoSpngjwrAjVlaY8GgkAgKYuaOaTma6+DXykUI07lX+joWnbe7PNngTh4JSgvRe2u+0DSN43orlpbL496vkmRn/mMyp/qtAIDzI4mDtXCYxm/2dGNDvQp7/+VqXN+sjaljEfIL5vvDzBnH5N7zvDXz/Y4xWFx+fGFr/ayPYzoa2Jq7SzTLlk/cvKIM/mAY701LwIw5fAiEaFZtEBi0ciGMCywJpze5Z1xXVRIBHN5gRtvM/cEwJrxB+INhMlOdAV4/PYwiPhfbmrWsvu5inQxOXxCDVk+sU8XhDWJ8HsmAmNG6MvmZO21UHX0kTzq1pjNi84KmQSp3KTBX5W4cwHEAJyb9OR79M3vj8SWIQswviMqdweyCRiacYnAdj5iReY5bZ471WSDgcbCicnYzTgatXIhAiIY1j/vHM0X7iANjDh+ubIrvdaCUZL+63DrWhTsAACAASURBVDpgAwC0j0wkDLhHJ7zwB8O4oUWH8pIiVCrEGHP4yNxdHsJ8f1RMcMeosOa5l9fOMyNQSwVzihfURDdFBhaCO28gBL3ZhSZtfrZkMqypVkAnF2FH28iU24dsUW+pDFTuNHLRgpq5o+mI4mjtdH/RaFUkk+fdycfesI0kxdgkEApj15kRXNesnXPPlCpLooqZj7zVgdYBG25fXQEA6B5L31rGYHaDooAqZfLHbKWiCAIeJzYzmm/EbBBI5S5p5gruegFcRdN03aQ/9TRN1wGYvcfqEqS0QNoB9WY3apNQUtLGvO5ye7E4qrdgdVUJhLzk2pqYdS+0eY7pHNNbZlwE9naNAwCuShDcqaTZGe6fTGu/DXwuBYc3mLAax1RJmO8lcxIfti3sz7AQYTapjAUJo8Ka6wr/bLj9QbzfMYaPLNOBO4c8eLGYD4WYjz4WRJl6xp0I00BTEmbiuYTDoXDj8jLs7RyHw3sxKTaYAQNzBo1MCLPLD38wt23/bDHu9MHtD824tjJJkExWtie/Ngnu2GX/BROs7gA+ynJLJnBRZGlH2wiWVxTj//vIYgCYV5DVb3GjTC5Ker8ERCw7GkqlGbFlYIO+6P6ABHfJM1dw9ziARMMJv2B5LQVPtVKCAasboTzveTeYXTGz3tlg/FhGc1iqn/AGcHbInvS8HTDZIHfhBgYuXxCfe+bYjOHrPZ1jWFomhyaBQb1aGqncZWsuY2zCi2G7F9cv0wGIVO/i0T+tVZjZTA6R1sy8w5ygcpfPM3cfdIzDEwglbQ9Qo5KwUrljLCKyZWA+H25aUQZ/KDxlNnYoGihkauYOyO/vTSroo8mA6ZU7VSz5kbmfc3Ib3xAJ7ljl9dPDKC7iY2uChOl8kIn4sWPrh7c0o7xYBKmQN8/KnSslGwSGfLC0cfmC6B6buYZDPSbo5CKilpkCswZ3NE3/hqbp0wnu+5/MLKlwqVWJEQjReZ05c/uDME745lTKBAClWAABl5NTO4Q/HdQjTAPbmnVJP+eiQe7C2DTE480zI3D6gmgbtOPMYETlzuEN4ITBiqsWJ74IKSUChGnA5slOy+qpaEvmXeuqQFGRttF46M1u8LkUyqMnbyZDN2i9dC0t8hVLtPLLzNyJBTyIBdy8nrljWjIZn8y5qFNLZtghhMN0ynO8naNO8LnUjA1/PrKmugQVJUV4c1Jr5pDVgxIxHxIh+5a4Cy0Jd7H7YOq19aLgULYqdwvj95kPeAMhvHNuFNtbdKz6PE7mjrWV+PymOqyrVYKiKCwqlaB7fH6Vu1TEVBiatDKM2L2w53Cc5Yn3LuDGX++fojofDIWx/4IJVzaV5q0HXz4yl6DKdyf9+85p9z083zenKOoGiqI6KYrqpijqe/N9vVzDZEvY9EhiG0YBLpnMDodDQSMX5qwt0+4J4A/7enHdUg2WJzlvBwCaBbZpiMeLR/tRoxJDxOfgL0cNAIAD3WYEw3TCeTvgYhbZkqXWzNaBSEvmulolalWShJU7g9mFKoU41jKnk4vA5VBEVCUPMbv8EHA5kE7a8EeMmvMzmZJKSyZDjUqMYbtnir3DPU8ewU/eSM2e44LRgXq1NK5ybb5BURRuXK7Dvi4T7J6I6fLb50aTauFPByYJl8+iKueG7Xh6f19SQb3eNNMGAbh4zs3k8cFUBcuLRXmdXC40juutcPlDuKEl+eRyqnx7WxN+eEtz7P+LNNK0K3dOXxAmpz+pzqzpMPN/XXEqZ9lib+c4/MEw3jl/sXvg9KANE94grpwlaU2YyVxXnLsm/fv70+67YT5vTFEUF8BvAGwH0AzgUxRFNc/+rPymVs2uylomiLWOJHnB1slFORNUeWp/Hya8wVk9qeIh5HFRIuYvWK+7LqMDJ/ttuHdDDW5dWY7XWofh8Aawt2scUiEPa2sSy7wzrXTZmo9q7bdhaZkcIj4XS8tkaE/gw2QwT8028rgc6OQi0mKUh1hdfiglgilZ1IhRc35W7lJtyQQi50eavlg5Hnf4cKjXHFPuTZZOoyPv5+0mc9OKcvhDYXzumaO47+mjUIgF+PnHlmfkvZgk3Hienqd9wRC++sJJ/MeO83j5xOCcj9ebXaiaZoMAABIBFyI+J6PHh9npg4jPQYNWRs6ZLHJ6MNJ5sroq8TWVbRo0UhgnfJjwpl5Bi403pKCUycCcp3I1dzfm8KIz2hY6uXtgb+c4OBSwaZE6J+sqVOYK7qgE/473/1S5HEA3TdO9NE37AbwI4KPzfM2copWJIOBxWPNHygTMHEmymR1dsSgnmVWry4+n9/dhe4sOy8qTr9oxaGWiBduW+eLRAfC5FG5fXYG719fA7Q/hn63D2Ns5hk0NqlmrBDFlwywEd6EwjbZBG1ZVlQAAmsvkMJjdUwQbgIjKXLw5gUpFEWnLzEMs0eBuMuosGTWnQ6otmcDF2U8mGXag2wQg4iGZbGumKypxvjjPlTIns7KyGJWKIpzst+H+zXV44+ubsbQsvq3KfFFJhOByqLyt3D29Xw99NOn0H2+cx8Ac13W9yR133IGiKKgkmfW6Mzv9UEmEqCgpIpU7FjkzaEetSsyacXkyNJRGzhfpiKqk43HHUF4sgkzIQ1eOgrtDPWYAwKYGFQ50m2KtmXu7xrG6WpHVz2AhMFcjPZ3g3/H+nyoVAAYm/X8QwPrpD6Io6gsAvgAAVVVV2LNnD5YuXQq9Xg+Px4O1a9fixIkT0Gg0EAgEGBwcREtLC7q6uhAKhbB8+XK0trairCyStR0ZGcGqVatw5swZcLlcNDU14ezZs6isrITf78fY2FjsNYuKilBbW4v29nbU1tbC4XDAbDbH7pdKpSgvL0dXVxcWLVoEs9mMLRIj+kdl2LNnD0pKSqBSqdDT04OmpiYMDw/D6XTGnq9SqSCTyaDX67P2Mw102VAjBk4dOZDUz6Txj2LUMYSJiTU4efJk1n6ml/edARXg4HoVFwcPHkz5c1pC92PIVo49e/bE/ZxsNlvs+fn4OSX67rWsXIXzJw/h1koFfBNmWLs7sV5D45V39kPjd2FzzdZZv3s1i5dhFW8I/d1F0Mt9Gf2Z+sbsaAz3Y1lpHfbs2YNiZxAKyo033tqNjauWxn6m2sUtaAz3Q+mhMDysiH1OtdQ4hsdMcDiWx/2cPjjZCZt9Avfdek3efU6znSMK9bvH3A+TEeViBfbs2RP7mXSuXnQ6pHl33mtevgpjXSexvVIN0/hY0p+Te8SIVbwh9AzXgDfWjkNdE9BwPCgP23Gu2wCv3TTnz9Q3akUR/BCZOnH+fKBgvntPfGwJLpxvw7KaEEaHBjL6ObWIHbB2n4JjY3leHU8dPXo8e8CB2ysCuLapCD854MUjz7yCr920DhRFzfiZGhsbIbF0okZVh/Pnz8/4nCpFfoSHzkGvL8nIz+Q0DWM5xwQtXwaNdwC733sf6y9fR8578/yZXPrTWFm7NKv7CJe+F3LKiyOHD8FpkKX0M3V2DGAVbwhKQQh79uxJ+XNaX+yApfsUHI6arH9OB4YEuLzIiNvLi3Cmx4MX/rELG1cuhn1Ej+2LiuBwOC6p714yP9Os0DSd8A+AEIAJAA4Awei/mf8HZnvuXH8A3AngyUn/vxfA/8z2nLVr19L5zv3PHqOv/9XeXC8jIXf/8RB922/2J/34P+7roWse2kHbXP4MrmoqJoeXXvLvu+hv/PVk2q/xnZda6Q0P72ZxVfnB661DdM1DO+h9XWOx2144bKBrHtpB1zy0gx6yumd9fiAYomse2kH/6p3OTC+V/uuRyLp6xhw0TdP0kNVN1zy0g/7Twb4pjzuuN9M1D+2g3283Trn90Xc66brv7aB9gVDc1//E7w7Sa/7jHTocDmdk/YT4bP3F+/TX/zL12PzFW+10/fffpEOh/PosjvVFvlvvnhtN+bkrfvw2/W//aKPD4TB92c/epdf97F265qEd9MFuU1LP/9vRfrrmoR1037gz5fe+VLj1fz6k733qSK6XMYNvvniKbvy3nbTB5KJp+uJn+eSHvXEfb5zw0DUP7aCfPdAX9/7PPXOUvvGJfZlaLr398X305585Sr9yYmDKOZeQPmMTXrrmoR30H/f1ZPV9A8EQ3fivO+mHd55P+bnff7WNXvWTt9N+7++/2kav+PHbWb+mhsNh+oqHd9Nf+vNxOhwO05sfeY++76kj9D9PDdI1D+2gW/utWV1PoQDgOJ0gXppLLZNL07ScpmkZTdO86L+Z/8+3RjoIoGrS/ysBDM/zNXNOrUoMg8WVNan5VIm0jiTfj63LgZH5i8cG4AmE8PVrGtJ+Da1ciDGHL28/h3SgaRrPHzagUlE0pf/81lXlkAp5aNJKY2qTieBxOVCI+VkxMm8dsKG4iI+6qFJgWbEIxUV8nJ+mmMm0vk1vJalUFCFMJ7bi6DO5YHb5MebIz7auhYrFObMtUyURIhSmYc+SCmuy9JkibUqLNKm3RtaqxDCY3eg0OjDu8OHu9dVTXnMuOkYdEPE5qFKm3iJ1qVAqyz8j8+N6C/5xaghf3FofG1+487JKXLdUg0d2deCts6MznpPoHMagkggy2gpvcvqglgpj53+imDl/2qLzdisqS7L6vjwuB7VqcXptmWY3quchgLRYK4PdE8j6NVVvdmPY7sXGBjUoisJNy8txoNuE11uHoRDz0VKR+mjOpU4uJbyOAWikKKqOoigBIuItr+dwPaxQo5bAGwjn5YbTFwxh2O5BdQqbDV3Ui2gkS4qZwVAYLxw2YFODCg2a9IUItHIRQmE6b0Ue0uGZA3oc6bPg/s114ExS/ZMKeXj0Eyvxo1uWJfU6SokgK0bmrQM2rKwqiQlvUBSF5jI5zk9TzDSYXeBQQKViWnAX3agM2mbOurh8wdgxdn44vkgLgX18wRAcvmBMmIchNsuZJRXWZNGbXeByqLTMb2tUEujNLnzYFZm3+8RlVRDyOOhNUqb8uMGC5RXFSSt0XoowSbh8gaZpPLyzHTq5CF++alHsdoqi8Oidq7CsQo6v/uUkXpkmsMKIqNUlsLxQSYUwu3wpW2kkQzhMw+LyQyUVxJQ6ydzd/GkbtINDAS0VmZk5nY1FpekpZhosLtTMI5nEmKpnW1Rlf3SmeXNDJGl90/IyBMM03usYw5bGUnIOTYOcBXc0TQcBfA3A2wDaAbxE0/S5XK2HLZiB6nxUzByweEDTF1U9k4Exms2WrcB7HWMYtntx74baeb3ORZnthZHBPKa34OGd7djWrMVnN9bOuP8jy3TY1JCcmpRKKsy4oIrLF0SX0RETU2FYWiZH5+gEQpMqqgaLG+UlRTN8hJhgL54dwuTqyfRgkZA5rK5IZU4pnS6owhg151cyRW9yo0pRlJYVQa1KjCGrB+91GNGoiVTF69SSpCp3dncAZ4fs2EgU3mZFKxfB4vLDFwzN/eAssKdrHCf7bfjGtY0QC6ZKEhSL+Xj+/vXYUK/Ed14+jecO6mP3JbJBYFBLBQiEaEx4g6yv2e4JIBimoZYKoZWLQFHEyJwN2gZtaNTIZnwPskGDRop+i3uKFctcBEJhDNu8aYmpMCxm7BCyHNwd7DahvFgU2z+3VMhjRYjZrJ0Iicmp+Q5N0ztpmm6iaXoRTdP/mcu1sAUjQWvIw+COWVMyHncMTHA3as9OdvXPhwwoLxbhuqWaeb3ORZnt/MkKp8uYw4uvvnASlYoiPPqJlfM28lRJBBmvaJ4etCFMA6urpwd3MngD4SnJD705fquwrji6UZkluONzKZwbtrO8ekIimMpcwspdngV3vSZX2gbiNSoJwjRwuNeCLY2RDUZ9aXLB3eE+M8I0sHFR8gqdlyLaPDpP0zSNX73ThSplEe68rDLuYyRCHp76zDpsa9biR6+fw5HeiMKfweyOa4PAoM6g113smJQKIOBxoJEJSeVuntA0jbZBe0r+umzSoJEiTKdWJBiyehAK0yl1Zk1HKRGgVCbMauUuFKZxsMeMTdGWTCBSKb9lZRm4HApbmkiCLB3y31m1wCgvEYHHoaDPQyNzZk2pzNwJeByopQKMTmT+YtE95sT+bhPuXl+d8CKZLNmuOGaSb/2tFRPeAH5371rIRfOXA86G4fRJgxUAsGaaPxAjqz7ZzDxigzDzgiTgRbzu4lXu9NEN9qYGNWnLzCLMrKZSIpxy+8XKXe436Qx01GIjUavcXEzucGA2GHVqCfotbgRC4Vmfe6jHDBGfg9XV2fPHKkSY8/SAJffByDvnjTgzZMc3rmmctdIr4nPxxF2rUFFShB+9fg7BUBh9pvjnMIaLbcvsJz/GHZHXLI0eg+UlRRjO0hjFQmXI5oHZ5cfKHAV3i6J2CKm0ZhoszNxn+jN3QMTMvMuYveDu/PAE7J4ANjdODeK+fk0jXv/aplgXFiE1SHDHMjxuZIC+P4ngbsTuwasn5zZHZQuD2QWZiAdFin4hWrkooagFmzx/2AA+l8In11XP+7WYC12+eigli90TwIFuM7505SIs0bHT+6+UCGHzBBCcY4M6H07229Cgkc7wpmnUSlHE52JP5ziASPuazR1IuDFK5HXXZ3KhrFiEtdUK6M1uOH3stzsRZnIxuJtauVOIBaCozFQm0mXM4YPbH0o/uItukgRcDtbXKQEAdWopgmE6bsJhMgd7TFhXq5zRakyYypoaBfhcCu93GHO6jnCYxmPvdqFeLcHtqyvmfLxYwMMPbm5Gx6gDzx0ywGB2zZo0LZUx1yP2r6MXK3eTgjsiqDIvzgxGukGyLabCsKhUCooCesaSr9xdiAZk8bwWU6FJK8OFMceU0YlMsu9CZC9wxbQuBxGfm5bHMSECufJkgBqVOKly+sM7O/Dtl05nrT/eEG1/S7Wtr6xYhNEMB0kuXxCvnBjETcvLYhfC+SDgcaCSCGB0FPZFjsmgrWTxIqOWCkDTgNWdGWVDmqZxqt+KNdUz1yzkcfHJdVV4rXUII3YPDJbZW4UrSoriHh+9pkhFprk8EvB2kLm7rJAouONyKCjFApjySMCIaZ9MpVNhMkqJADIhD5fVKmJzN0ygOJuoypjDiy6jk8zbJYFcxMeWxlLsPDOaEbGRZNl5dgQdow48eF1j0l0jH1mmxZVNpfjvtzvg8odm3VSXySOzeJlIkpocF9sygYvnzFz+Pgud04N28LkUlpSlL+o2H4oEXFSUFKE7SfEmAPjwggn1agk08vlVuhbrIqMT/ZbMd5/Z3H48+WEvNtQrSYWOZUhwlwFqVRIYzO5ZT67jDh/eOjsCALG+/UyTqP1tLrRyUcbbGz+8YILDF2SlasegkeefzHaqML3vTTr2LjKqaEtdpuwQ+kwuWN0BrEnQknb/5jqE6Yj651ytwpUKMUbs3hlVRr15anBHRFWyg8XlB4cCSopmVv+z0e6bCkxwl27ljqIo/Oz2FvzLRxbHbquPvtZsc3eHeiLn800NZN4uGba36DBk8+DMUO5mZ587qEd9qQQ3ryhP+jkUReHHty5DOHpqmm22U17EQxGfi5EMBHfm6DGpEEeCu/JiEfzB8IJSis42bYM2LNHJIeRxc7aGBo00rrDJgMWNN9tGptzmDYRwuNeMKxfPX3xkcVQxsz0L19T/frsTE95g0krfhOQhwV0GqFaK4fQFZz25vnR8AIEQDRGfgyO9loyvqcvoQL/FHZO6TQVdVNEsFeWmVDnVb4WAy8GaGvYqVFq5sODbMrtGHZAJeSgvZi+rxVRdMrURPxGdt1tbEz+4q1KKcfOKMvzlSD/ORjd0iYbAKxRFCIXpKT6LVpcfNncAdWoJdHIRFGI+mbvLEmaXHwqxYIoVB4NKknkV1lTQm1wQcDlzej/OxkdXVUyZm1NIBFCI+eidI7iTiXikpShJtjVrweNQ2Hlmpn9csoTDNO558jBeOGJI+bnGCS+OG6z46MqKlCXX69QSfPHKelAU0DjLtZWiKJSVZGa8weT0QSkRxNZeTuwQ5kU4TOPMoB0rcjRvx7ChXoVOo2NGIulnb57HV/9ycsroz5E+C3zBMCvKkot1MqilQjy8sx1jGex8ahu04S9H+/GZK2pjs/gE9iDBXQZgBvETKWaGwjT+cqQfV9SrsLmhFEf6Ml+5+9mb7ZAIefj0hpqUn8sYmY9lMFA62W/Fsgp2M2VamSijJ6ds0DnqQJNONm+FzMmoMzjcD0Tm7eQiXmwoPB5f2FoPpy+I5w7qoZULUSSI/7kz/mSTFTN7J1VkKIpCc/lM7zxCZohnYM6glgnzSlClz+RCtUrMukdSnVqCvvHEwd2BHhM21KuIN1OSlIgF2Nigxq6zI2m3Eh7ps+BAtxmvtw6n/Ny3z42CpoEbl+vSeu9vXdeEt7+5NaENAkNZsSgjQicmpz8maASQ4G6+9JldcPiCrI5CpMNtqyrAoTBFl2Hc4cN77WMAgL9Pun1P5xiEPA421M+/W0DE5+Kpz1wGs9OPzz1zLCPz7OEwjR+8dg4qiRDf3NbI+usTSHCXEZj5Ib0pfs/yns4xDNk8uPeKGmyoV0JvdmdUsOSDzjHs6xrHg9c2JtyYzQYT3I1mqMXRHwyjbdCesI0vXTRyIcYdvqwNBrMNTdPoNDpi3jNskenK3al+K1ZVK+JWdxiWlRdja1MpfMHwrOpe8bzu9NPa7ZrL5OgYdWRUIIYQweJKHNypJIK8qtz1mWYXuUiXOrU0YVvmgMWNAYsHm4gFQkpsb9HBYHannaR5+fgAgIjx9FxKptPZdWYUDRrprJW32eBwqKQ6YsqKizJWuVNN8p1kgswhIqqSFm+fi1SQ19crc7oOXbEImxtL8erJIYSje5hXTw4iGKbRoJHilRODsdv3do1jQ70KIj47yfGVVSX4v3vWoGPUgS8/fwL+IHvXVpqm8du9PTg9YMO/3bSEFQVwwkxIcJcBKhVF4FCJK3fPHzZAIxNiW7M2lmnJVPUuEArjZzvOo04twX1X1Kb1GrrogO5IhuSV20cm4AuGMxDciRCm80vBLxWMEz7YPQEsYTm4KxELwKHmrtz5g+GU2x0d3gA6jY64YirT+dLWegCzq3uVRRMLk0VV+kwucDkUqqKtnM3lcviDYfTMUk0hsIPZNXUjORm1VACHL5jR9u1kCYdpGCxu1JeyH9zVl0owOuGFK05Gm5m329hAxFRS4fpmLThUJNBKFYc3gJ1nR6CTi+AJhNAxkryMu8npw5E+M25sSa9qlwplxZHZdbaTUOZplbsSMR9FfG5Klbu2QRt+v7eH1XUVIjRN4+8nBnF5rXLelgJscMeaCgzZPDjcawZN0/jb8QFcVqPA169piN0+YHGjd9zFutn31Us0+Pnty/HhBRO+8ddTrAR4/WY37nnyCP777U5sa9bitlVzK9MS0oMEdxlAyOOivKQo5jsymQGLG3u6xnHXuirwuRwsLZNDJuLhcIbm7l44bEDPuAv/duPStGW5mcpdpkRVTvZHPdFYnLcDAK2ssO0QOqNKmenMSc4Gl0NBIZ7byPy5g3rc+OsP8c65xBsuuzuAl44NwO2PbHRPD9hB04nn7SZzxSIVvnhlPe5YE98wGIi0iGhkwil2CH1mF6oURTEvquayyGzE+RFiZp5pZq3cSTMr1JMKw3YP/MFwhip38UVVwmEa/2wdgloqRKMmcUsyYSYqqRAb6lXYeSb11swdbSPwBsL44S3NAC5eT5LhnXNGhGlg+/KylN4zHXTFkWTjOMvJRpPTFxPJAiLzfeUlopSCu8O9Zvx8V0deHLu55GS/Db3jLnx8beJrUjb5yDIdZEIe/n5iEMcNVvSOu/DJdVWR20U8vHxiEHu6IlYCV7EgpjKdT6yrwg9vbsZb50bxxT8fn1fi7vXTw7j+8b1oG7Tj4duX4/efXsvquAlhKiS4yxC1KgkuGJ0zLlT/t6cbFIC7Lo+oQnI5FNbVKjNSufMGQnj8vQvY1KDCtUs1ab+OTMSHRMDFqD0zQdLJfhvKikUoK05f+CAehW5k3jkaqZotZjm4A5JTNnz9dGR+5buvtM1oJ6JpGq+1DuHaX+3Bd19pw1dfOIlAKIwTBisoClhVNXegTlEUvr99KdbPMSfQoJHiSJ8l1l7bNz7VmHpRqQQCHoeIqmSYUJiGzROAUpy4LRNAXrRmMi3xk43I2SJRcPfMQT0O9pjx4LUNZNOSBtuXl6HX5MLHfnsQH//tQdzz5OFZLScYXj4+gEaNFNtbdNDJRTFBp2TYdXYEdWoJ690R8SiPXt/YVMx0+4Nw+0NQy6YekxUK8ayKrtNpiYr/nM2hYmk+8PcTAyjic3HjiswH+8kg4nNx88oy7Do7iqf390Eq5OGmFWUQ8bm4ZWU5dp0dwZttw6hSFqWtCjwXn99ch4dvX449XeO4/7ljsURuKjh9QfzwtbNYrJPj3W9vxd3rq2cd2yDMHxLcZYitTWqcH5nAi8cGYre932HEX48O4IEt9VMU3NbXKdE77mJdtv9Qjxk2dwAPbKmf92ZDWyzC6ERm2jJPGqyst2QCk4K7AhVV6Rx1QiMTQpHGnORcKCWCWbO0BrMLZ4bsuHt9NfzBML71t9ZYcHV+eAL3PX0UD77YioqSInz9mgZ80DmO7796Bif6rWjSyCBjsY/+3g01MJjdeOtsxAtLb3ZNkR3ncTlYopMVpKhKIBROeUYoV9jcftD0TI87BqZylw+iKn2mSFCQiQ1PvODu7JAdj+zqwLZmbVqiVQTg1hXluL5ZC7GACyGfg0M9Zvzj1NCsz+kec+Jkvw13XlYJiqKwpqYk6cqd1eXHwR4ztrfoshKMMx0wIyzOwjGJFLVkqjfsuhoFOkYdSR+LjLLr2eFLN7jz+EPYcXoENy4vg1TIy/VyYtyxphKeQAi7zo7ilpXlMc/NO9dWwhsI43CvBVc2lWb0O3z3+mo8eudKHOox44n3LqT8/BcOG2BzB/DjW5pZT+IT4kOCuwxx/+Z6bGlU40evncOZQTtMTh+++/c2LNHJ8J3rm6Y89uLcHbut2EAG+QAAIABJREFUmbvbjZAIuNjIwnB/WbEIfabZvfvSYWzCiyGbB6uTmNFKFbVUAIqKVHreazfipzvO493zRtbfJ1N0GidYF1NhUElnl61/80zER+crVy3Cj29dhkO9Zvz49XN44LljuPHXH6K134af3LoMr35lE75z/WI8eG0j/n5iEPu6xllvr71+mQ71agl+u7cbxgkf3P5QzG+MoaWiGKcH7PAFcz/vlQqfe+YYvvlia66XkRQxA3OpMO79pXkV3LlRxOdCmwFjXBE/YjD8QecYDnabYHH58Y0XT0Eh4eORO1aQql2aFIv5+MN9l+GFBzbghQc2YFVVCT68YJr1OS+fGACXQ+H21ZE2ujXVCgxaPUl1a7x73ohQmMb2luxUaS5W7thLkjLH2vTKHeN3tn+O3x9DsZiPKmURzg0VXoKMLd4+NwqHL5g3LZkMa2sUsbn0u9ZVxW5fVVWCRdGZ4qua0u/MSpaPranEpgY1dqe4h/L4Q/jjh73Y0qieYitDyCwkuMsQXA6FJ+5aDbVUgC+/cALffuk0JrxBPH7Xqhly/8vK5ZAKeay2ZtI0jd3tRmxtKmXFXmDbUi3aRybwfscYC6u7CJNlzcRBz+NyoJYK8eT+Ptz/3HE8tb8Pj77Tyfr7ZIJQmMYFozMjLZkAUBmdCX2tNX5m/M22EayuLkGlQow711bilpXl+PNhA47prfj2tibsf+gafGZjbUzu/ZvXNeJTl0cuPGtr2FUZ43IofPHKepwdmsCfD+sBRBQLJ7NtqRZOXxAHuzNvK8IWZqcPB3pMeLfdGFecI99gZjRVCSt3mbXYSAW92YUalThjrT83tOhwesCGu588gjU/fRd9Jhce++SqtNSICfHZ3FiKtkEb7O5A3PuDoTBePTmEqxdrUBqdr14TnfU9mURr5ptnRlCpKEJLRXY8tjJhZG5yMsfk1IRLS3kxlBIB9kbnsZJheUXxgq/chcI0usfiC+68fGIAVcoirK/LrUrmdCiKwoPXNeJjqyumeO9RFIXPbqyFUiLAFVlS5716sQY9464pHntz8dej/TA5/fj6NcTyIJuQ4C6DKCUC/OaeNTBOeLGvaxwP3bAES3QzLyQ8LgdraxSsiqqcHZqAccKH65ZqWXm9ezbUoL5Ugv98s51VWdyT/TYIuJyMXWC/dV0TvnTlIrzwwHp89epF6DQ6YPfE3yzkEwazC75gOGOVuy9ftQhraxR48MVWPL67a0pFts/kwrnhCdwUFRmgKAqP3LEcv/7Uahz43jX4xrWNKBZPbbukKAo//WgL/u+eNbg5A/MKt62ugFYuxB/29QKYOUu1sUEFqZAXk7EuBPZ0joOmI6qkc1Uo8oFY5S5BACMWcCHic/JCnVZvcmVsBgUAfnBzM1p/dD2e+ew6fPHKejzysRXYuIgoZLLJlkY1wjRwqDf+sbG3axzjDh/uvOxipWVZuRwCHmfO1kyLy4/93SbcvKI8a5XWTBiZm2OVu6nBHYdDYUujGh9eGI/J5c/FsvJiGMzugrg+pstj73bhul/tw99PDE65vXfciYM9ZtyxpjIvZ8FuX12JX31y1Yzv6qc31ODIv14LSZbaSK9eEqkQ7ulKLsnvC4bw+309WF+nxOV5FjQvdEhwl2FWVyvwq0+swuc31eFzG2sTPm5DvQrdY07WTvzvthvBoS4ejPOFz+XgBzc1o9fkwvOHDay8JhDJsLJtXj6Zu9dX43vbl2BTgxqbGtSgaeCEITPKpGzSORrJLmYquCsRC/D8/etxx5pKPL77Ar75t9aYEtbOaEvmjZMU5MQCHm5dWT7rLAKPy8GNy8tY89qZjJDHxQOb6xEI0RDwOLEWp8n3X7NEg3eirVaFwPudYyiVCSEX8fBee/63C89VuaMoCipJxFsylwRDYfRb3FPmMjOBXMTH1Us0+P72pfjEpHYpAjusqiqBVMhLmPh4+fgg1FIBrpl0jRPyuFheUTynqMpbZ0cRCtMZSUTNRlmxKCNtmfGOySubSmFy+pOeRW6piFSFzi3Q6p3dHcCzB/Xgcyl875U2HOiOfK+6xxy458kjkAp4+MRlhXUcUxQVU43OBnVqCWpVYnyQRAdXOEzjmQN6GCd8pGqXA0hwlwVuWVmOH97SPGtGaFtzpMK26+wIK+/5XrsRa2sUrLYJXbW4FFubSvH47i5Y02y9CoVp9Iw74fYHI+blQ+yblydidZUCfC6Fo33Jq6nlik6jAxQFNGoyp+Im4HHwyztX4F8+shivnx7GHb89iEGrGzvaRrC2RjFF9Ccf+NT6ashFPNQmaLe7oUUHi8uPY/r8D94DoTD2dY7jmsUaXLVYg/c7xvI+KLVEW8BmE/ipL5Wg0zi3wmEmMVjcCIbpjFbuCJmHz+VgQ70S+7tnBndmpw+72424bVXFjM3t2hoFzg5NzDp/+8bpYdSrJVhWnp2WTAadvIj1tkyZkBc3obalMTJ3l2xrJvO7WKhzd88e1MPpC+KFBzagvlSCL/35BF49OYg7f3cIgRCNF7+4Ie+uefnIVYs1ONhjTmiLMGBx45G3OrD5kffxX7s6sL5OiU0N2WkbJVyEBHd5QoNGiiU6Gd5sm39wN2zz4NzwBGstmQwUReHfb1oKlz+En+9qn9Hu4Q2EcKTXjNdah/CHfT34n/cu4L12I8YdPlhdfvxubw+2/uIDXPvoXjT/8G1s/K/34c+AeXkiigRctFQUF8Tmv3PUgRqlGEWCzFQ0GSiKwlevbsBTn7kM/RY3bnziQ7SPXGzJzCekQh4e++QqfPcjS+Lef2VTKYQ8Dt46m/+tmcf1Vjh8QVy9RIPrmrUwu/xoHcjvpIPF5YNcxJs1U7y6qgSdoxNpyWWzQShM4ydvnAefS+GyJLwWCfnN5gY1DGb3jBmff7YOIximcWecSsua6hL4Q2GcTRCkjE14cbjPjJtXZq8lk6G8RIQxh481I3OT0xebdZ1OqUyIZeXypIM7tVSIsmLRgpy7c/qCePpAH65bqsXldUo887nLIRJw8e2XTkMi5OHvX7oiphhKmJ2rl2jgC4ZxqHfmfPuQzYPbfnMAf9jXiyadDE/ctQrPff5yIjKVA/JH75WAm1eU4ZfvdGHY5plXBolp8bqumd3gDogYat+/uQ5/2NeL0QkffnnnCmhkIrzfYcSPXj+HAUv8lhMOBYRpYEO9El+9ugEWlw+94y5MeIPY3Ji9WZXLa5V4+kAfvIFQRtoH2aLT6MhYS2Y8rlmixetf24wv/Ok43H4Xti/XZe29U+HaWRIWEiEPWxpL8fa5Ufzolua8vqB80DkGPpfC5kY1QmEaPA6F3e1jrIvRsIl5FgNzhpVVJQjTwJlB+5z+hfMlFKZxTG/B2hpFLOD81bud2Nc1jodvX476UmIkXuhsaYpUnz7sHsc9qojFBE3TePn4AFZWFsc9RzLJwlP9VqyNE+BHjNKBW3LgZaYrFiEUpjHu9LEiCW92+qFOoF4LAFubSvHHfb1weANJ2dMsKy9ekF53zx82wO4J4GvXNAAAKkqK8NznLsezB/vwnesXx2yTCHOzvk4JEZ+DPR1juHrxxZZobyCELz9/Ar5gGLse3IKmDInBEZKDBHd5xE0ryvHLd7qw88wIHthSn/br7G4fQ51agkUZ2tx8f/sSVCvF+P/bu/O4quvs8eOvwyoiogIuyKKAW26471tZWTaVVlZmi9lo+zLTOjVTs9RU35n21TEtK23SzEwbtVxyj9xQQRQBFUQFRNmU5XLfvz+48EMFvChwL7fzfDx8gPdz7+dzrn7gc8/n/X6f8/el8Vzz1np6hfizZl8mUa2b8dGUvkS1bkab5k1wEyH+aC6xqafILijmhuj2DZqwVGVAh1Z8vC6ZnamnKlpQOJuM3EIOZhXwu17BDXrcjoG+LHl4OEdzzjTaXjTjerTlp73H2ZWWQ287Gqk7yuqEDAZ1DKhYwziwYyt+ij/OM+OqHpV0Btl2JHflzetj007Ve3K3eMcR/rgglrBWTXnsik74eLnz/pokbhsQyuRBYfV6bNUwIgJ9CfZvwobELO4YVJbc7TmSS8KxPP5+Y48qX9O6eRMignyZH3OYKYPDz7uJt3TXUbq29aOTAz58Vm5kfqm/Y+PTc9mRepLxPau/TozqHMSHa5PYlHSCq7tf+IZdz/b+rEooq97bUEU66lthSSmzbKX4oytdEy4Lbs7rN/d2YGSNUxNPd4ZFBrJmXyYvGYOIYIzhz4v3sCsth5l39tPEzgnotEwn0tG2BmDpJUzN3Hs0l81JJxjbrf76nogIUwaHs/SR4QT5ebM5+QTPjOvKD4+OYFyPdkTZmlj7enswoEMr7hsRwdPjujo8sQPo36HsTu6vddxTsCYJx3J5emEsszekkJp94RLCn2xMAWBi3/b1Hdp5fLzcG/WIx9hurXF3E5Y7cdXM1OzTHMjIP6vY0RXd2pCYkc+hEwU1vNKxypK76kcJoKx/YmgrH3amnqr3eNYlZtKiqSd+TTz444JYHvxyO71DW/DXG7rX+7FVwxApG93elHSiYk3q11tT8fZw4/re1Sc1f72+O0mZBby+/OzWN+mnzrD10El+V8Nr61NdNTLPLihm+udbaeHjxTPXdKn2eX3DWtLM28PuqZk92jfHGOwuwtIY/PfXVC3FX8dGd23N4ezTJGcVkJFXyDurDrBgWxqPXB7FVXbcRFD1zzVuzbiQ63oF89ryBFKzTxPaqumFX1DJ9sMnuWd2DK18vbhrSIf6CbCSTm38WPrIcAqKS/H3ufCUD2fQoqkXXdr4EdNA6+4WbkvjhcW7sRootqTxt6Xx9GjfnPcn9yU84PyCDzlnSvhyy2HG9wqucruqWYumXgyJCGDJznTuHxXplOdlea/IylX+xnZrzd+XxvPT3gymDe/oqNBqlF1QTO+QC4+GRoe2ZFs9/3wZY9iUdIJRnYN469ZoVsQdZ/meozxzTdd6q7yrHGN4pyC+3prGE//dycnTxcSkZHN197Y1/myP6BTE3UPCbeusWjM0qmzqf/ma9oauklmuLhqZW0qtPDxvOxl5RSyYMYTWftVPKfTycGNIZADr9mdibKMsNSmvmLnnSA4DOjjvFHF7WUqtzNqQTL/wllqKvw6Ntk2XnvTR5ooqymO7tebxsZ0dGZaqREfunEx5IYvycvT22nQgiymzfqGlrxcL7h9S68TwYnm4uznlB+iaDOjYku2HTtbZovaqFFusPLNwF08uiCU6tAUbnhnD2idH8/y13UjJLODV/yVU+bovthwiv8jC/aMuflrub92DoyPJyCvkrtkx5BY6R8+mHYdP8p91yfzt+3g+2ZBCx0Dfs6o5hgf40ql1M1Y4aTGYklIrJ08X06qa4g2V9Q7xJz2nkIzcuqsKeK7EjHwy84oYFhmIiDCuR1veuq1Po51OrKo3IiqQ5k08WLMvg9wzJYzv2Y6nrq5+tKrcs9d0IyLIlycXxLJmXwYPzdvOa8sTiA5t4bAbZ+WNzC+l5dErPySwKekE/5zQ066p549e3on3Jve1a9+t/bwJbOZdbTGaxmZF3HFSs88wfaReT+tSaKumTOjTnsuCm/PMuK4seXgYM+/sj7sT9gj8rdKROycTFtCU3iH+fL8rnXuHd6woFJCZV8QvKSdIP3WGWweEnZVQbT98kns+/ZWOAb58Pm0grXVxcI0GdGjFF1sOk3Asr+JOZbmkzHxyz5TQ55wKnst2HWXZ7nTemBRtVyGW2RtT+O/WVB4eE8XjYzvh4e4GfvD7kRHkFVl4Z1Uie47knHX8wpJS5mxMYVTnIK3cdQmGRgXy/uS+PPjldu6ZHcPcaYNq7M9X3wpLSrnzkxjyiyz4ernT1r8Jv69iTe3EviG8tjyB/cfznG7NwqakE5SUmrPWrFSnT1jZc3aknrJrnc/F2GDrfTZUS2y7vJa+Xuz4y1W4CbUqkuTj5c6bk6KZ+OEmps75FX8fT6YO68C04Y77oC8itl53F5fcfbMtjdkbU5g6rAM39Qu58AuAniH2X0tEhB7tm7tErztjDDPXJdEx0LfOK4crePPWaEeHoGqgyZ0Tuq5XMC//sJfOL/yPAF9vfLzczqpC+WP8cebeOwgfL3cycgu5//NttG3ehPnTB9dpXztXVT49IyYluyK5sloNszemVKzR+P6R4RVrBI/nFvLsN7vIK7LQISCRpy9Q9CK/yMLHPycxqnMQT1Zxh/m+ER35bNNB/rVyH59OHVjx+IJtaWTlF/PA6Mg6eZ+/ZVd1b8t7k/vw0Lwd3DvnV+ZOG+iw6qir9maQX2RhztQBjO4cVO0H1FsHhPLmT/v5fPOhaotFOMrS2HT8vD0YZZuOU5Puwf54uAmx9ZjcbUrKIjygKSEtG2aGgnKsix0R6B3agvcn9yW/yMJ1vdo5RYXkdi0urpF5bOopnvt2N0MiAvjTtd3qIbIyPdv7sz4xi/wii0NvilVn2a6jHMst5K4h4We1ZTmRX0R2QXFFoZyYlGxi03L4x409dERJ/ebotEwnNGVwOK/f1ItHL+/E2G6t6dW+Bc9e05XFDw3j7dui2XroJA/P287pYgsPfLmdvEILM+/qp4mdndr5+xDS0ofZG1P42/fxfL01lbvnxPCPZXsZ0SmQ5j4ePP7fnRVNcF9aEkdxqZXLu7bm43XJ7E6r+a7mnA0pnDxdwh+urHr+efMmntw/KpK1+zLZalubdPjEaT5am0SfsBYM0rUBdWJcj3a8eWs0MQezeWHxHoxxTJPw72PTCfLzZmSn6hM7gFa+XlzXqx2LtqeR5yTTSQGKLKUsjzvGVd3b2vXhuImnO93aNa+3oiqWUiu/JGczNLLhWqioxmtcj7bc3C/EKRI7uLhG5pl5Rcz4fBtBzbx5/46+NfaavFRjuram1GpYuDW13o5xsWatT+ahedv5+9J4bnhvI3uO5FBQZOHtnxIZ+foarnprHX9fGs+Z4lL+sz6ZVr5e3GznCKdSrsT5bssofLzcmTTg/OasUFZqPLfQwp8X72Hsv38mPaeQ9yb3oWvb5g0cZeP2+NjOfL7lEPNiDlFYYsXH051XJvTk9oGhrNqbwX1zt/Lmj4n0C2/J//Yc46mruzBlcDhXvvEzTy2MZcnDw/HyOP8Cm3O6hJnrkxnbrU2N6yHuHhrOJxtSeH3FPgZ3bMVH65LxcBP+dUtvp+7P1thc3zuYpIx83l6VSLd2zRu8WEluYQmr92UweWCYXXeP7xrSgUXbj7Bo+xHuHtqh/gO0w7r9WeQVWvhdb/uLUPQO9WfxjnRKrabO75rvOpJDXpGFYTolUzVClRuZe9iRpFmthofnbefUmWK+eWBovd/E7RvWkr5hLZi98SB3DungFKNexhjeXpXIWz8lcm3PtlzXK5gXl8Rxw/sb8ffxJLugmGt6tKWVrxefbEhhZfwxUrPP8PjYTk6T1CvVkDS5a4TuHBxOdn4xb/60nxmjIriugfuhuYKb+4Vwc78QSq2GlKwC/H08CfIrK/M+9rI23D4wlI/XJdHSVl1z+sgIPN3deGVCT+6bu5V/rdzHpP4h+Hh54O/jWTF9ZdaGZPIKLdWO2pVr6uXBw2Mieen7eGJSsrkhOpjnrulWUSpb1Z3HruhEwrFcXl4WT0SgL1Gtm5GZX4Svl8dFtecwxvDrwZN0D25+wV5QK+OOU2yxcn20fT+j0aEt6B3iz+dbDnHXkHCnSPS/j02nZVNPhkXZP1IWHdqSL7YcJikzv87XD246YFtvpyN3qhGq3Mi8qMTKgm2pDI0MrPbna9GOI/ySks1rN/VssLXYvx8RwQNfbmdl3DGu6emYyqKVvfHjft5dfYBb+oXwz4k98XB3Y1hkIK+vSOBYTiEPjomqaFg/vlc7nl64i6Ze7tw5ONzBkSvlGJrcNVKPXhHF+F5t661R+W+Fu5sQ1fr8f8MXxl/GxgMnSD15mll396+YBjP2sjbcGB3MzHXJzFyXXPH8iEBfokNbsCLuGNf2bMtlwRceSb19UBhZ+cWM6BRY7w2ff8vc3IQ3JkUz8YNNTP3017O2fTSlL+N61O7Dy4q4Y9z/xXa6tPFj1t39a6xMuyQ2ndBWPvSpRUP1O4d04MkFsfy8P5PCklK+2X4EAd6+rQ8+Xg17F/p0sYUf449zY5/2tZoKVl54ZWfqqTpP7jYcyOKyds11GrpqlMrbITy3aDcbErOwWA3vr0liXPe2PD++21m/T/IKS3j1fwn0CWvBLf2qns1TH67q3pawVk2ZuT7Z4cndhsQs3l19gEn9Q3h1Yi/cbCOJ/k09eXlCz/OePzQykB+fGMWpM8UENKu5L6dSrkqTu0ZKRIhq7VwV9VyJr7cHn08bSHJWAX3PqZz5f7f05oY+7ckrtHCm2EJGbhGxaTmss12on7Cz14u3h3uVBVdU3Sv///x+11H8vD0I9PPirZ8SeXrhLnq097e7MEd+kYWXlsTTMdCXozlnuP69DXw4pR+Dq0jOT+QXsfFAFjNGRtRqBO66Xu14eVk898wpS0QDm3mTXVDEg19uY+Zd/et1vc25VidkcKaktFZTMqHsZodfEw+2HsxmUv+6+1B6priU7YdOcfdQvSOvGqfgFmXJ3YbELCYPCmP6yAi+25nOe6sPsGZfBk9d3YV7h3XEzU14d/UBsvKL+OTu/hVJTUNwdxOmDe/Ii0vi2HYom37hjlkHnnO6hKcWxhIZ5Mtfr+9h97+Bj5c7Pl7aFkX9djkkuRORW4CXgG7AQGPMVkfEoVRNwgN8q+yH5Onuxpgurc973BhDcalVmyg7qdbNm5y15i4yqBnj39nAo/N38N8ZQ85Lmtbtz+STDSnMGBVRMQXw3yv3cTyvkEUPDMXfx5P75m5lyqxfuKV/CBP6hNA/vGXFB5Afdh+l1GrsnpJZromnO8+Pv4xNSVlc3zuY4VGBLNiWxnOLdvPkgljenBTdYB/0lsYeJcjPm0Edazey7OYmXNmtDV9vTaN7sH+drR/cnJxFcam1oim1Uo1N5zbNePPW3vQJbUkHW6/Lh8ZEMbFve/7yXRz/WLaXjQeyeGhMFHM2pjCpf4hd/ezq2i39Q3jjx/38Z10K/e50THL3lyV7yMwrYtGDQxt81oJSjZmjRu72ABOBjx10fKXqnIhoYteIhAf48s+JPXlk/g7+b8U+/nhVZ7w93MkrLOGVH/YyPyYVDzdhfWImT4/ryrDIQD7bdJA7BoVV9EFc/NAw/rE0nsU70pkfk0pISx+6tPGjWRMPth8+Sec2zS6q2FH5mtBytw8M4+TpYl5fvg+L1TCue1u6tWtOh4CmdhVluBg5p2tXDOZcr0zsSV6RhReXxJFXWMJDY6IuaQ3hsZxCnlu0m2D/JgyuZbKplLMQESb0Ob+CYzt/H2be2Y+5mw/x8rK9rNmXiZ+3B09dXXPrnfrS1MuDOwaF8eHPSWxOOsGQyPr/mTt0ooD9x/MptVo5kJHPdzvT+cOVnekV0vDJrVKNmTiqPDiAiKwFnrR35K5///5m61Yd5FNK1Z3nFu1ifkwqItDGrwklpVZOni7m9yMjmD4igr98F8ey3Ufx8XTH19uDVX8chb+P51n7KCiysDL+GMt2HeVoTiEFRRYKikv545WduW1gWJ3EaYzhjR/38+HaJCzWst/bAb5ePD2uC7f0C63z0byZ65J45YcEfnh0hF1rSKtSUmrl6YW7+HbHEW6IDubWAaEM6hhQZbK4/fBJkjMLqixdfrrYwqSPN5OSWcDCB4bSrZ1WB1auKy49hxcW72HywDBuqcNpzbWVc6aEmz7cxPHcQr55YGit1s8eOlHAs9/sJiOvENuvK4ZGBjB1WMcq17lvTjrB1E9jKCyxVjzWP7wlX00fXG83sJRqzERkmzGmf5XbNLlTSv2WFVusLI87RkpmAYezT5NfVMKMUZEVay2NMczeeJB/rdjHvyf15loHFxgospSSlFHA3qO5fPXrYX49eJLeIf785Xfd6RvWotrRsW2HTjJrfTL3jYioqCwHUFhSyvrELEZ0CqwoG15qNYx8fQ3tW/rw9YwhlxSv1Wp4bUUCczcd4kxJKYHNvJnQJ5gZoyIJtBU8+GLLIV5aEofFavjyvkFnVQ60Wg0PfLmNlfHHmXVXf67o1uaS4lFK2S/t5GkmfLAJL3c3vn1wKK2bn1/RuchSipe7W8Xvnt1pOUz9NAaL1TAsKhA3EYpKSlm7P5Nii5VRnYO4d3hHRnYKRESIScnm7tkxhLby4dWbetHEwx13NyEiyLdB1xgr1Zg4JLkTkZ+AtlVset4Y853tOWu5QHInItOB6QBhYWH9Dh06VA/RKqVUzeqjZ9ulMsaweOcRXvkhgcy8IsJaNWVstzaMvaw1Azq0wtPdDWMM82IO89KSOEpKDW4CM0ZF8tgVnfhp73Fe/V8CaSfPMGVwGP+4saz63Iq4Y8z4fBsf3tG3zqrlnS62sCYhk+9j01kZf4wmnu5MHdaB7IIS5sccZkyXIJIyC/DycOOHR0dU9JF8bXkCH65N4oXx3bhvRESdxKKUst+eIzlM+ngzHQJ8+eCOvhVrBYstVv69ch8z1ycTEejL+F7BhLdqyp+/20PLpl7MnTbwrIreJ/KLmPfLYeZuOURmXhFRrZtxfe9gPvo5iXb+Tfhq+pCKlkRKqZrpyJ1SSrmwvMISvtuZzqq9x9mYdIJiixV/H0/GdAnCAN/tTGdU5yBentCD91Yf4KtfU/Hz9iCvyELXtn5EBjVj2e6jfD5tICM6BTH5P1s4mFXAuqfH1MuUqKTMfN76KZHvY9MBeGhMJH+4sgtr92Uw7bOt/OnarkwfGcmCrak8tXAXkweF8fKNPZyi759Sv0Vr92XwwBfbKSm1MnlQGL/rHcxLS+KIS8/l+t7BZOYV8UvKCawGurVrzmdTB1Q5ygdlSeGy3enM3nCQ3UdyiAj05avpg6t9vlLqfJrcKaXUb0RBkYX1iVn8GH+c1QnHOXm6hIfHRPHElZ0rRh5gvh7sAAAJdElEQVTXJGTw4c9J3NS3PTf3C6Wk1Mr4d9ZzuriUd2/vw80fbeaZcV15YHRkvcaacCyXkwUlZxVruO+zX9mcdIJXJvbkyQWxDOoYwJypA3R6llIOlpFXyDurEpkfk0qp1dCyqSev3dSLq7qXTdLKzCtiZ+ophkQG0Mz7wvX6jDHsOZJLSEsfWmrfSqVqxemSOxGZALwLBAGngJ3GmKsv9DpN7pRSyn6lVkPumRK7PjjtTD3FxA824uPpjsVq2PLcFQ75wHX4xGnGvvkzxRYrkUG+LHpw2HkFbJRSjpOcmc/3sUe5bWAobXS0TSmHqCm5c8itUGPMt8aYEGOMtzGmjT2JnVJKqdpxdxO7E7To0BY8ODqKguJSJvRp77A76WEBTXnqqi60b+HD7HsGaGKnlJOJCGrGY2M7aWKnlJNyVJ87pZRSTubRKzrh18SDCX3aOzSO34+M4L4RHXWNnVJKKVVLmtwppZQCwMvDjRmj6nednb00sVNKKaVqT1eoK6WUUkoppZQL0OROKaWUUkoppVyAJndKKaWUUkop5QIc2ueutkQkEzjk6DiqEAhkOToI5dL0HFP1Tc8xVZ/0/FL1Tc8xVd+c6RwLN8YEVbWhUSV3zkpEtlbXa0KpuqDnmKpveo6p+qTnl6pveo6p+tZYzjGdlqmUUkoppZRSLkCTO6WUUkoppZRyAZrc1Y2Zjg5AuTw9x1R903NM1Sc9v1R903NM1bdGcY7pmjullFJKKaWUcgE6cqeUUkoppZRSLkCTu0sgIrNFJENE9jg6FuV6RCRURNaIyF4RiRORxxwdk3ItItJERGJEJNZ2jv3V0TEp1yQi7iKyQ0SWOjoW5XpE5KCI7BaRnSKy1dHxKNcjIi1EZKGIJNg+lw1xdEzV0WmZl0BERgL5wFxjTA9Hx6Nci4i0A9oZY7aLiB+wDbjRGBPv4NCUixARAXyNMfki4glsAB4zxmxxcGjKxYjIH4D+QHNjzHWOjke5FhE5CPQ3xjhLDzLlYkTkM2C9MWaWiHgBTY0xpxwdV1V05O4SGGPWAdmOjkO5JmPMUWPMdtv3ecBeoL1jo1KuxJTJt/3V0/ZH7/ipOiUiIcB4YJajY1FKqdoSkebASOATAGNMsbMmdqDJnVKNgoh0APoAvzg2EuVqbNPldgIZwI/GGD3HVF17C3gasDo6EOWyDLBSRLaJyHRHB6NcTgSQCcyxTS+fJSK+jg6qOprcKeXkRKQZ8A3wuDEm19HxKNdijCk1xkQDIcBAEdEp5qrOiMh1QIYxZpujY1EubZgxpi9wDfCQbdmMUnXFA+gLfGiM6QMUAM86NqTqaXKnlBOzrYP6BvjSGLPI0fEo12WbYrIWGOfgUJRrGQZcb1sT9RVwuYh84diQlKsxxqTbvmYA3wIDHRuRcjFpQFqlmS0LKUv2nJImd0o5KVuxi0+AvcaYNxwdj3I9IhIkIi1s3/sAY4EEx0alXIkx5jljTIgxpgNwG7DaGDPFwWEpFyIivraiY9imyl0FaBVzVWeMMceAVBHpYnvoCsBpi9t5ODqAxkxE5gOjgUARSQNeNMZ84tiolAsZBtwJ7LatiQL4kzHmBwfGpFxLO+AzEXGn7Gbf18YYLVWvlGpM2gDflt0PxQOYZ4xZ7tiQlAt6BPjSVikzGZjq4Hiqpa0QlFJKKaWUUsoF6LRMpZRSSimllHIBmtwppZRSSimllAvQ5E4ppZRSSimlXIAmd0oppZRSSinlAjS5U0oppZRSSql6JiKzRSRDRC7YrkNE3hSRnbY/+0XklD3H0OROKaWUUxORT0XkZjufa/eFsw7i2mT72kFEJld6vL+IvFMPxzsoIoF1vV+llFIN5lNgnD1PNMY8YYyJNsZEA+8Ci+x5nSZ3SimlGj1brz6oxYXzUhljhtq+7QBMrvT4VmPMow0Rg1JKqcbDGLMOyK78mIhEishyEdkmIutFpGsVL70dmG/PMTS5U0op5TRE5C4R2SUisSLyeaVNI0Vkk4gkl4/iichoEVkjIvOA3VD1hbOKY3wqIh/ZLqL7ReQ62+NNRGSOiOwWkR0iMsb2eHcRibFNjdklIp1sj+fbdvkqMMK2/QlbXEttz2klIottr9siIr1sj79kG2Vca3tP9iaDT9liiRGRqJrej1JKqUZhJvCIMaYf8CTwQeWNIhIOdARW27MzjzoPTymllLoIItIdeB4YZozJEpFWlTa3A4YDXYElwELb4wOBHsaYlFoergMwCogE1tgSpYcAjDE9bXdOV4pIZ+B+4G1jzJci4gW4n7OvZ4EnjTHlSeLoStv+CuwwxtwoIpcDc4Fo27auwBjAD9gnIh8aY0ouEHeuMWagiNwFvAWUJ3LnvR9jTKH9/xxKKaUamog0A4YCC0Sk/GHvc552G7DQGFNqzz41uVNKKeUsLqfsApYFYIypPAK32BhjBeJFpE2lx2MuIrED+Nq2v0QRSaYs0RpO2boGjDEJInII6AxsBp4XkRBgkTEmsRbHGQ7cZNvnahEJEBF/27ZlxpgioEhEMoA2QNoF9je/0tc3L/B+dtYiTqWUUg3PDThlW1dXnduw3Xy0d4dKKaWUMxDAVLOt6JznlSu4yGOdexxzzn7//wZj5gHXA2eAFbYROHtVtc/yY1d+T6XYd8PV2PF9VX9XSinlZIwxuUCKiNwCIGV6l28XkS5AS8puMtpFkzullFLOYhUwSUQCoGy9Wj0e6xYRcRORSCAC2AesA+6wHbszEEbZdMkIINkY8w5lU0J7nbOvPMqmVlal8j5HA1m2i3m1RGSViLSvZvOtlb5WvthX9X6UUko5ERGZT9nv7i4ikiYi0yi7RkwTkVggDrih0ktuB74yxth9w06nZSqllHIKxpg4EXkZ+FlESoEdwD212YftwjkaCBSRNOBFY8wnVTx1H/AzZVMh7zfGFIrIB8BHIrIbsAD3GGOKRORWYIqIlADHgL+ds69dgMV2Yf7UFne5l4A5IrILOA3cfYH43YAoqi8K4y0iv1B2c/b2mt5PTcdRSinV8Iwxt1ezqcoqz8aYl2p7DKlFIqiUUko1eiLyKbDUGLPwQs9taCLSA7jXGPOHWrzmU5z0/SillGpYOnKnlFJKOQljzB7A7sROKaWUqkxH7pRSSimllFLKBWhBFaWUUkoppZRyAZrcKaWUUkoppZQL0OROKaWUUkoppVyAJndKKaWUUkop5QI0uVNKKaWUUkopF6DJnVJKKaWUUkq5gP8HLjLeMBYWT84AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot eigenvectors to confirm successful eigenvector decomposition.\n",
    "\n",
    "plt.figure(\n",
    "    figsize=(15,2)\n",
    ")\n",
    "\n",
    "loc_eig = bioframe.slice_bedframe(cis_eigs[1], 'chr1:10M-60M')\n",
    "plt.plot(\n",
    "    loc_eig['start'],\n",
    "    loc_eig['E1']\n",
    ")\n",
    "plt.axhline(0,ls='--',lw=0.5,color='gray')\n",
    "plt.ylabel('E1')\n",
    "plt.xlabel('chr1 position, bp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:57:03.649044Z",
     "start_time": "2018-09-25T16:57:03.513680Z"
    }
   },
   "outputs": [],
   "source": [
    "# Digitize eigenvectors, i.e. group genomic bins into \n",
    "# equisized groups according to their eigenvector rank.\n",
    "\n",
    "Q_LO = 0.025 # ignore 2.5% of genomic bins with the lowest E1 values\n",
    "Q_HI = 0.975 # ignore 2.5% of genomic bins with the highest E1 values\n",
    "N_GROUPS = 38 # divide remaining 95% of the genome into 38 equisized groups, 2.5% each\n",
    "q_edges = np.linspace(Q_LO, Q_HI, N_GROUPS+1)\n",
    "\n",
    "# Filter track used for grouping genomic bins based on bins filtered out in Hi-C balancing weights\n",
    "# Doesn't do anything with eigenvectors from the same Hi-C data (hence commented out here),\n",
    "# but important for external data, such as ChIP-seq tracks\n",
    "#eig = cooltools.saddle.mask_bad_bins((cis_eigs[1], 'E1'), (c.bins()[:], 'weight'))\n",
    "\n",
    "# Calculate the lower and the upper values of E1 in each of 38 groups.\n",
    "group_E1_bounds = cooltools.saddle.quantile(eig['E1'], q_edges)\n",
    "\n",
    "# Assign the group to each genomic bin according to its E1, i.e. \"digitize\" E1.\n",
    "digitized, hist = cooltools.saddle.digitize_track(\n",
    "    group_E1_bounds,\n",
    "    track=(eig, 'E1'),\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:57:03.982931Z",
     "start_time": "2018-09-25T16:57:03.651273Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 0, 'chr1 position, bp')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3UAAADQCAYAAACgN0jgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9eXhcd33v//7Ovs9IM9LYWr3JCYm3OCZxFjuQpJSyBAiUQpM0baHQ3tLS9f6g/fVy6aVPeUop99db+tymUBpiytYSoNASICWxsziJs1h2EkcjOx5t1mzS7DNntu/vjzPfo5E0+5wzmiN9X88zj6TxLF/Jc875frb3m1BKweFwOBwOh8PhcDgcdaLZ6AVwOBwOh8PhcDgcDqd9eFDH4XA4HA6Hw+FwOCqGB3UcDofD4XA4HA6Ho2J4UMfhcDgcDofD4XA4KoYHdRwOh8PhcDgcDoejYnhQx+FwOBwOh8PhcDgqRqf0GxBCtADOAJinlL6DENIP4JsAdgC4DOD9lNLleq/h8Xjojh07FF4ph8PhcDgcDofD4fQmzz//fJhSOlDt3xQP6gB8HMCrABzlnz8B4FFK6WcJIZ8o//z/1HuBHTt24MyZM8quksPhcDgcDofD4XB6FEKIv9a/Kdp+SQgZAfB2AF+quPtdAB4sf/8ggHcruQYOh8PhcDgcDofD2cwoPVP3vwH8dwClivu8lNIrAFD+OljtiYSQjxBCzhBCzoRCIYWXyeFwOBwOh8PhcDjqRLGgjhDyDgBBSunz7TyfUvoApfQIpfTIwEDV1lEOh8PhcDgcDofD2fIoOVN3C4C7CCFvA2AC4CCEnAAQIIRsp5ReIYRsBxBUcA0cDofD4XA4HA6Hs6lRLKijlH4SwCcBgBDyJgB/RCm9lxDyOQD3A/hs+ev3lFqDkvz9Y9PQEILfvG33Ri+Fw+FwOCri9XAKD5y8iEKR1n2c1ajDf3/rVbAYOrtUf/O5GewbduLaIWdHr1PJ+fkYXphZxq/ctEO21+RwOJxWSAoFfO5HF5DOFdt6voYQfPDGMRwada26/68feQ0TXhvedWhYjmV2jW6oX67lswC+RQj5EIAZAL+4AWvomCd8YeQKJR7UcTgcDqclHnzqMr753Cy2OUw1H1MoUQQTAo7u6sdb921v+70KxRL+9OHzeOfBIXzhlw61/Tpr+fS/v4znLi/jPdcNw27Sy/a6HA6H0yw/Or+IB5/2w+swQktIy8+PZwv40cuL+N5v34IdHisA4NtnZvF3P5vG/TeN86CuGpTSxwA8Vv4+AuCObryvkthNOrweTm30MjgcDoejMk76Qjg2MYAHf/2Gmo8RCkXs+9QjODsX6yioW4hmUShR+IKJtl9jLRcW43jusmgve24+hpt3e2R7bQ6Hw2mWU74QPDYDnv7EHdBoWg/qZiJp3PXFJ/Dhr57Bd/7bzfAFEvjTh8/jlj1u/Nk7rlFgxcqitPrlpsVu0iORLWz0MjgcDoejIuaW07gUSuHYRP1AyKjT4qptdkzORTt6P/+SmHycDiZRKtVv92yWE6f9MGjF7cPkXEyW1+RwOJxWKJUonvCFceseT1sBHQCMuS34+3sO43I4hd868Tw++tDzGHKZ8MVfPgydVn0hkvpW3CPYTTrVB3V/8K2XcOQzP5Fuf/dfvo1eEofD4WxqnvCFAQC37W2s6nxgxIXJuVhHwZg/kgYAZPMlzEczbb8OIykU8PAL83jHwe0Y6TN3HHRyOJzuEU4KeNcXn8SlUHKjl9Ixr1yJI5LK4XgT59J63Lzbg/9517V4cjoCIV/Cl+4/ApfFINMquwsP6trEYdIjKRRQlCnzuRE89loIbqsRP3/tNmg1BE9Mhzd6SRwOh7OpOekLYZvDhD2DtoaPPTjiRCJbgH8p3fb7zVQ8V44WzIdfnEcqV8R9R8dxsBx0cjgcdfDUxQjOzkZxdhMkY076RA/rW/d03v5979FxfPbu/XjwQzdgz6C949fbKHhQ1yZ2kziOmBTUWa0rlSii6Rx+7hov/uI9+3FwxIVoOr/Ry+JwOJxNS7HcLnRswgPSxFD/gRFRka2Tapg/koLXYQQATAU6y85TSvG1035cO+TAoVEXDow4MbecQSQpdPS6HA6nO0zOiueSzbDfOzUVxtXb7BisIzjVCh+4YQyHx/pkea2Nggd1beIoq30lsuo8MBLZAkoUcFnE36PPYsByOrfBq+JwOJzusZzKIV8sde39JueiiGcLTbcLTQzaYNJrcHa2/WqYP5LG/mEnBu1G+DoM6s74l3FhMYH7jo6DELISdM7zah2HowbYsRrLqHPvykgJBZzxLzXVxr6V4EFdm7BKnVrn6lgA11fuG3ZZ9VhO50GpettJORwOp1mKJYrbP/8Yvvq0v2vveXIqDEKabxfSaTW4dsjZdqWOUoqZpTTG+q3Y67VjusP2y4dfnIfNqMNdh4YAAPuGHSAEmOwg6ORwON2hWKI4Xw7q1F6pe+b1CPJFimMTPKirhAd1bWKXKnUqD+qsK5W6XKGETL49A0cOh8NRE0upHJbTeVzsomDAKV8IB4ad6LM2P4R/YMSJlxfiKLRRUQwlBaRzRYy7LdgzaIMvmOwocffiTBSHx/skM3S7SY9dHisXS+FwVMDFUFIy6VZ7pe7kVBgmvQZHdqi7XVJueFDXJiuVOnUeGCyoYwo/feU2zKUUb8HkcDibn3B5Diyc6M48WDybx4uz0ZYzywdHXMjki5huI/icKStfjrktmPDakM4V21bAzOSKmAokcHDEuW59k/Mx3uXB4fQ4Z8vzdDajTv1BnS+EG3e6YdJrN3opPQUP6tpE9e2XKfGAltovy1/VXpLncDicZgiVg7lQl0Q+npqOoFiiLctvHygHUe20ODI7g/F+CybKim6+YHuVyVeuxFAsUewfXh3UHRhxIpQQsBjPtvW6HI6a+MHkAmIq3Sedm4/BatBi/7ATURVrKDTr9bkV4UFdm9hVLpTCKnX95WCuv9wOxMVSOBzOVkCq1HUpqJuci0KnIbhuzNXS83a4rbAbdW1JkPsjKWgIMNJnwUTZQmG6TbEUJtZycHT1+g+Uf+5EzIXDUQOBeBYf+5cX8e3nZzd6KW1xdi6GfcNO9FsNqq7U/WDyCgDgzVcPbvBKeg8e1LUJq9TFVVqpi6bz0JCV34O1Xy6rNAPF4XA4rSBV6hJCV1oHF+NZeB0m6LWtXXY1GoL9I862/OD8S2kMucww6DTosxrgsRnb9qqbnIvC6zDCu0Y+/JrtDug0hM/VcTY9gXI1+kpMfVXpXKGEVxfiODjqgsOsV21QVypRfO0ZP27Y2Y/dA429PrcaPKhrE5NeC4NWg7iKK3UuiwEajeiVtNJ+ySt1HA5n88MqdNl8Camc8gJRwbiAwbJfXKscGHHhwmIcQqG1dfojaYy7LdLPE4O2tr3qJudjkoVBJSa9Fnu9dm5Cztn0sHNGQIWtxq8tJpArlnBgxAmXRY+oStXOH/eFMLuUwb1Hxzd6KT0JD+o6wG7SqXemLp2TPOoAwGXmQikcDmfrEKoQSAl1QSwlEM/Ca2/PJPfgiBP5IsWrV1qrsjE7A8aE14bpNhQw49k8LoVSOLBmnk5a36hou1AsrX7dV6/E8VsnnkemC0Ezh6M04YS4PwrGu9OyLSeT82Il/cCwCy6zHoUSlZQw1cTXTvvhsRnw1mu3bfRSehIe1HWAqoO6VF4SSQFEPyS7SceFUjgczpYgnMxVfK/8Jm0xnsU2Z3tB3fXjomz36UuRpp+TyOaxlMqtq9QlhULLoibny1W4A6PV5wGPTwwgni3gsdeCq+7/u/+axn+eX8ST0+GW3o/D6UWYqJIaRYEmZ2Pos+gx2m+Gs5zEj6qsBXNuOY1HLwTxgTeOwaDj4Us1+F+lA+wmvaqFUiqDOkAUS+FCKRwOZysQTgoY7TcDUL5Sl84VkMgW2m6/HHSYcPU2O05OhZp+TqXyJWPCW1bAbLEF8ywL6mpU6u68xotBuxEPnV4xcg/Gs3jk5UUAoj8fh6N22HkiEM+qrnXx7FwU+0dcIIRIXVpqU/H8+rMzIAA+eOPYRi+lZ+FBXQc4zOqt1EXTeUkcheGyGLhQyiZCbRcdDqebhBIC3rDNAUD5Sl2g3K7VbvslABzfO4Azl5eRzjV3zZGCOndF+2VZAXMq0Fob5+RcFKP95pqm6XqtBh+4YQyPT4Ukb7xvPDeLQoniKq8dp3xbq1LHz72bE1apEwolxDPq2ftlckX4gknJY9IhVerUk8QXCkV887lZ3H61F8Mu80Yvp2fhQV0H2I0qr9StuUD3WfRcKGWT8Knvncf9X3luo5fB4fQkhWIJS+kc9nrt0BDlK3VMWKHd9ksAODbhQa5YwjOvLzX1eP9SCoBoPM5w24zotxpwsUUj88m56iIplXzwhlFoCMHXnvWjUCzh68/O4NiEB7/0xlFcCqcwu5Ru6T3Vii+QwNV/9iO8Hk5t9FI4MhOuOE+oqQXz5YXVHpMus7j3i6uo/fLUVBjhZA73HOVVunrwoK4D1DpTl8kVIRRKq4RSANGInAulbA4m52N49vUISiWeMeZw1rKUyoFSwOs0od9q7EKlTtwAettsvwSAN+7oh1GnaboFcyaShsdmgM2oW3X/sMvckiR7JClgPpqRsvy12O404843DOLbZ+bwn+cXcSWWxb1Hx3F8r2gQvFWqdefmYxAKJbwebk9llNO7hJICBuziMawmBUymTMs8JtneT00aCrPLYlLoUIPk0laHB3UdIM7UqS+oWypX49bO1DGZW476CcYFZPMlzEczG70UDqfnYG1UAzYDPDZD1yp1az3eWsGk1+LGXe6mgyN/JI2xink6htdhxGILQR3bEDaq1AHAfUd3YCmVw58+fA7bnSbccfUgdg/YMOQ0bZm5uoXyOVeNewNOfcIJAfuGxJZtNVXq1npMqlEoJZwUoNMQae2c6vCgrgPsJh2SQmGdjHOvs5xiQd3qg6PfYkBSKCBXKG3EsjgyUSpRaRPZrtEwh7OZYUHcgN2IAbsRoaSyHQqBuACLQbuuatYqxyc8mA4mpcCBUoonp8NV5+xmltKr5ukYgw4Tgi0EsWfnoiAE2FdDJKWSm3e7sdNjRTxbwAdvGINOqwEhBMcmBvDEdBiF4ua/trBEWpwHdZuKbL6IeLaAa4fE4yCoqqBudfu0xaCFXktUZUAeSghw21a8lTnV4UFdB9hN4gU6Kajr5M2qcesqdeUZOzUNz3LWs5TOoVBONLRrNMzhbGaYnYHHZsSAzbhqVkYJAvEsvA4TCOlsQ3J87wCAFTXJrzx5Gfd86Rl854X5VY8TCkUsxDKr7AwY2xwmLKVyTRuZTwUS2OG2NhWQajQEH7p1J2xGHT7wxlHp/mN7PUhkC5KK5mZmPipu9tU6b8+pTqScDB/pM8Nl0UviR71OLJPHpfBqj0lCxIqXmjqzwsmc1PrKqQ0P6jrAYRIrXWoaNgUg2RZUE0oB1NVnzVlPZa9/q9LlHM5WgFXqPDZWqRMUVSwUg7rONyQTgzZ4HUac9IVxciqEz/zwFQDApdBqUY7ZpQwoRdWgjq2jWQPlYFxoae333DiGZ//0DgxWtJreuscDQraGtQFvv9ychCvOGV67STXtly/PV/eYdJr1qtq7hhICPDYe1DWCB3UdwCp1ajt5M4XLakIpALhYisphQZ3dpMM0b7/kcNYRTortkFajDh6bEblCCQkFOy4CcaGjeToGa2U8ORXCx/7lBez12rHTY8XM0uqgjv081l+9/RIAgonmNqXhpICBFqwYCCGwGFZX9VwWAw6MuFry2VMjlFIpqFPThpnTmMqWba/TpJr2y1oek06zXlVdWeGkgAEe1DWEB3UdYC9X6tTWZrGUEtfLZG0ZK4pI6jnQOethbSE37XLDF0xyzyQOZw2VWV+P3SDdpwSUUizGs9gmQ1AHiC2YiWwBOq0G//grRzAxaMPlyGq7gBWPuurtlwCwGGvu9xX/VtX96VrhtgkPXpqNqmqOp1VimTzSObGtVW3JXk59mEKux26E125UTfvl5FwUY/2WdZ1ZLotBNV1ZlFKEkwI8vP2yITyo6wC1VuqW0znYjDoYdKv/+1mljhuQqxtWqbtljwfpXBELLSjdcThbgXCFNPmATQxylJqri2XyyBVKq9oRO+FNVw3g2IQHD9x3PUb7LRh3WzCzlF5lX+KPpGE1aOGuYhbOKobNSLKncwWkckVZZlkOjrpQosClFj3y1ESl2rDakr2c+qy0bBvgdZgQSgqqEMmbnIthfxU7EpdZr5oESyyTR75IeaWuCXhQ1wFSUCeo48BgRNM59FnXy8KuBHW8UqdmAvEsPDYD3rBdlF6eCvAWTA6nksrqk1SpU8irjmX05ZipA8RZ7oc+dCOO7OgHAIy5rcgVSghUtFP6IymMua1VhVn6LHoYtJpVj69FOLEiKNMprFKglupAOyyURVJsRnV62HJqE04KcJh0MOq08DpNKJYoIgr7W3ZKuI7HpMOsR0wlx6IUUPNKXUN4UNcBK+2X6jp5L6fz65QvAcBs0MKk12zqi+5WgM3vTAzaAADTXCyFw1nF6kqd+FWpSh0TVJCr/XItO8otlv6KFkz/Ulq6fy2EEAw6jAg0UcGX/Pxk2ExthaThfNkg+aptdtXtCzj1CVW0/3klA/LeDurO1fGYdFn0SAgFVdiMsPOQHG3gmx0e1HWAWtsvo+kcXFWCOkC88HKhFHWzGBPl0/usorEy96rjcFbIF0tYTuel6lOfxQCthihYqevceLwe42UxlJlyUFcsUcwtZTBWI6hja2lmQyqJQ8hRqSvPbG/m68tCLAuDToNxt4W3X24ywomcdBywY7lVBcw//vZZfP7Hr8m+tlrU85hkJt5q8FNk56FBXqlriGJBHSHERAh5lhBylhDyMiHk0+X7+wkhPyGE+Mpf+5Rag9KY9FoYdBrEVXbyFit169svATY8u3kvuluBYGJFPn1i0A5fkFfqOBxGpOxRx6pPGg2B22qQWg3lhlXElPJYGnKZoNMQ+MuKl4vxLHLFkhTsVWObw9Rc+2VyRca9UxwmPTRkc7dfzkczGHaZ4TDpVZfs5dSnUqhjm7P5uVRGUijg4Rfnu6oAOzkXw+4BW1WPSTUJ41X6inLqo2SlTgBwO6X0IIBDAN5KCDkK4BMAHqWUTgB4tPyzanGY1Nc7v5zKVW2/BMRsKhdKUS/5YgnhZE7KJE54bZgOcAVMDodRLVDx2IzKVeoSWfRZ9DDptYq8vk6rwXCfWVLA9EfE4K6a8iWj6fbLcobcLUPbk0ZD4LIYNnX75UI0gyGXCQ6TDslcYZV4DUfdhBIrkvpuqwEa0lpQd/piBIUSlczplYZSism5KA5UmacDVip1ahBLCSUE6LVEWjOnNooFdVSElQj05RsF8C4AD5bvfxDAu5VaQzewqywjly+Kfkw1gzqrui66c8tpZPPFjV5GzxBMMFGGclA3aENCKKjGKJXDUZpKvynGgN0oBXtyI5dHXT3G+i1S+yWbrRvrr99+mcoVkWzgzRdOCuiz6KHXyrNVcFn0m7pSt1Cu1NlNelAKJHPq2RtwapPNF5EQCtI5Q6fVYMBubCmoO+kTK3ThpNCVPcuVWBbhZA4Hq8zTAYCzbGkVVUFQF06KFjTVhJ84q1F0po4QoiWEvAQgCOAnlNJnAHgppVcAoPx1sMZzP0IIOUMIORMK9a5hqd2kU1XvPLugVlO/BMRKnVouusupHN7yhZP4m59MbfRSeobAGlGGCa8dAODjYikcDoAK8Y+1lTqFhFIC8aziQd0Ot1Wq0Pkjaei1BEMuc83Hr3jV1d+UVgrKyEHfJq7U5QolBBMChlxm1c7bc6pTbba02blUxilfGFqNGJQ0Ou7k4MnpMADg0Gj1oI61X8ZVENSFEvKehzYzigZ1lNIipfQQgBEANxBC9rXw3AcopUcopUcGBgaUW2SH2FXWfsn6p+sJpUTTOVW0jfzr83NI54p49NXARi+lZwiWg7pBaaZOVMDkc3UcjsiK31RFUGc3IJwUFGlTFoM6ZTck424L4tkCoukcZpZSGO2zSBvIarDzQ7BBpaHSpF0O+iz6TSuUEohnQSnKQR1Txu79DTOnMSvG4yv7pkG7qelK3exSGq+HU7jjarGGsVDhZ6gUX3tmBrsGrA3bL9WQxGeVOk5juqJ+SSmNAngMwFsBBAgh2wGg/DXYjTUohd2oV0Wmg8Hm5eoJpZQoel78pVSi+NozfmgIcDGUWmX6upVhGUBWGXDbjOi3GjDNFTA5HADiBsFm1MFsWJlxG7AZkS9S2edLCsUSQonutF8CYpXOH0nXVb4EKgzIG4ilhJM5WTPkohBXb19b2mVuWbwGDfNK3aZjpVK3chxvczbffslaLz94wxgAYE7h/cr5+Rhemo3i3hvHa7Ysqimoq/QV5dRHSfXLAUKIq/y9GcCdAC4A+D6A+8sPux/A95RaQzdQW6WOZUnrCaUA6CmxFEopHnzqMuaWV3yYnrwYxuVIGr95224AwBO+3m3R7SaB8kBxf8X/755BG2+/5HDKVNsgsMBF7rm6SCqHElXOzoAx7haVLi9HUpiJpDFeZ54OqJBkj9X/feWu1PWrbGa7FVj1ZXX7Ze9cRzntI6kvVlTqvHYTltN5CIXG83Enp0IYdplx8x43AOUrdV97xg+TXoP3Xj9S8zF6rQZWg7bnhVJKJYpISt7k0mZGyUrddgA/I4RMAngO4kzdDwB8FsDPEUJ8AH6u/LNqEYVSevugqGSl/bLWTF3vGcRejqTxqe+/jPv/6VmpgvjQ0364rQZ8/M4JeB1GnJwKb/Aqe4NALItBuwmaitar8X6LlEXmcLY61ebE2KxMUOa5OqU96hisUvfSbBQJoYAxd207AwCwGXWwGXV1Kw0poYBMvihzpU4PoVBCJrf5xK3YRn2701TRfqmehC+nNizZ47ZWzNSVbQ2CDebqCsUSnpqO4NiEB0adFgN2o6JBXTybx3dfXMC7Dg43VIt0WQyIZnpnr1eNaCaPYony9ssmUVL9cpJSeh2l9ACldB+l9M/L90copXdQSifKX5eUWkM3sJt0SOWKKKpgBg1YqcD1W2urXwK95V3iC4itgxdDKfzu11/E3HIaP301gPe/cRRGnRbHJgbwxHRYNf8HShJIZKV5GYbDrK7EA4ejJNWqTx6pUifveY+1Q29TOKgzG7TwOow45ROTWzsatF8CgNdhRLBO+2W12cNO6cWkoVwsxDLw2Iww6bVwmMVKnZpGMzi1CSUEuCx6GHQrW2aphblBC+bZOTHRcnyvqA0x7DJjQUFbg+88P4dMvoh7j443fKzT3PvjQ9XUijm16cpM3WaGtVkkVZKRi6ZzMOg0MNfwTGLtl0up3jnQmcjHJ3/hajz2WggfeOA0KIBfLvenH987gFgmj3PzsQ1cZW8QiAvrNpBqSzxwOEpSbU6MVerkVsAMSBYjym9IxvutmC6fK+t51DG8DlNdFT5WnZBX/ZK192++oG4+msWwSzz3OsqVurhK9gWc+oSTwirlS2DlmG6kgPn4VBgaAty8W2y9FIM6ZSp1lFKceGYGB0dd2F9DIKUSp7n31c6r+YpyarPeZp7TEisn7zycNVoae4nldA59Fn3N4VmmitlLlbrpYBLbnSZ89LbdmI9m8NWn/bj96kGMlluObt3jASFi33ot+d6tQiCWxa17PKvuY5/RZLagis8op/fJFUr4o2+fxZuvHsB7rqs9t9FrZHJFxDL5dRsEp1kPnYbIH9TFstAQUbBIacbcFjx7eQmEACN9zQV1z75eu1FmpVInn0ABu74s91DSUC4WohlJbdio00CvJbz9cpNQrbrPkqef+v55fP7Hr9V8biCexYERl/TZH3KZ8NNXA6CUyu679tzlZUwHk/jc+w409XiXRS8lgnoVXqlrDR7UdYjaVK6WUvmaIikA4DDpoNWQnsqk+oIJ7ClfLP/sHdfAbTXi7Qe2S//ebzVg35ATp3wh/O4dExu1zA0nJRSQEArr5nfYZ1QtiQdOb0MpxZ999zy+f3YBAFQV1P3n+SsAgOvH+1bdr9EQDNjl96qbW05ju9Nc115ALpg4yjaHCaYanRiVeB0mBBNZlEp01QwuQ4lKHWv776XrixxQSrEQzeC2cosdIUR18/ac2oSTAg6sMfF2mvX42Jv3wL+UrvEskWuHnXhfhWDJkMsMoVBCJJWTvfr06pU4AOC2q5qzAXOa9T1vPs4rda3Bg7oOcZjV5UcTTedqiqQA4sXIZdb3jPplqUQxHUzinhvF/nC9VoOP37k+cDu+14P/+/glJLJ5aUh9qxGs0eplr6gmczid8pUnL+ObZ2ZBiHrOe4wTp/3YNWCVWqEqEc2E5Z118S+lm2qFlANmYzDWQPmS4XWINg7L6VzVSmIoIYAQrFLS7RR27emlThA5iGXySOeKqwzf1aaMzalNtUodIQR/9PNXtfxa7DOyEM3IHqgsxrPQaQg81uZe12nRI5bOK1I1lItQQoBBq4HDxMOVZuAzdR2itkrdcjpXUySF0Wc19MxFd245g2y+JLW11OLYxACKJYqnLka6tLLeY61HHcOhss8op3c5ORXCZ374Ct5yjRc37OhX1Wfq5YUYXpiJ4p4a3k1eR/O+U83ij3QvqGO2Bs2+34rQQ/XqZCiZg9tqgE4r3zbBZWaVOnUlAxqx4lG3cu4Vgzrlf898sYR8saT4+2xVMrkiUjn5VGCHK4I6uQnEsxi0G6tW3qvhNOuRK5aQzffu5ydUVivu1aCz1+BBXYdI0sWCOi5Sy+m81Ntdiz6LXvKz22h8ZdPsCW/9oO7wWB/Mei2euaRqMdWOYEp269svubw2Rx4+8W+TmBi04wu/dAhOs15Vn6kTp2dg0mvwvsPV20XlrtQlsnkspXIY669vLyAXO91WaDVEalVvRCP1vnBSXo86ADDoNLAZdZuu/bLSo45hN3bn+Pi9b76EP/jWWcXfZ6syX/6/lTuom1dAATMQz2KwBaVdlmTpZVsDbjzeGjyo6xA1VepyhRKWUjkMNjg5uSyGnlFEYsqXewbsdR9n0GmwzWmS3TxYTax4Yq1tv+RGuJzOoZTiSjyLt1zrhdWoU9XMkOjdNO3LEsYAACAASURBVI+7Dg7VnCv1OkyIZwuyeaj5I+KsTTP2AnLgtOjx8H+7uSkpc6BSva/65jKUWO/nJwcuix7LPZI0lAumcrrNubZSp/y+4FIoJdn+cOTn6Uti98/aOdx2cVn0MOu1ClXq1qtfN1oLgJ42IK+mVsypTc0mVULIOQA1NdAppc3J62xyJBGKHj4oGKEkm7mqf9D3WfSYnOuNi64vkMSg3diUwIdDBUO/SrIYE2AxaGEzrj6s1ZR44PQu6VwRlALW8udLTTNDD78w39C7qbJytcPTeXVtpiygMNaloA7AOjGHegza67dfhpMCdsrwd1hLv9Ww6dovA7EstBqyypy6W0mPRDYPodC77XNq59RUCCN9ZtmSM4QQDLlMirVf3lJlXrgWzJy8V5L41QglBBxswp6BI1Jv8vAd5a+/Xf76UPnrPQDqy/1sIYw6LQw6jSo2N80a4fZZxItuLwzPTgcT2OutX6VjuLZ4UBdIZLHNYVr3f7bSftn9v81jrwVx/XhfXfGapy9GJNUuAHjTVQPYNdBcCxmneyQF8RzHgjqHSYdkrlBTPbFXoJTiodN+HBxx1g162HlxUaagjlXq2Kxbr2HQaeC2GrBYpVJHKVWs7UnsBOmNpKFcBOJZDNiMq1ROu5X0SGQLSAmFnrhebzbyxRKeuhjBXYeGZP3bDrnMUlsnAFwMJZHJFbFvuP3gJZ0rIJEttNR+2etBXbFEsZRSpmNgs1IzqKOU+gGAEHILpfSWin/6BCHkSQB/rvTi1ILDpFOFyWiwfPEebGCE22c1IFcoIZMvwmLYOMUhSil8wSTef2S0qcc7zXr4IymFV9W7zETSq2Y6GAadBsYNSDwsxrL41a88hz9529X4yPHdNR/3sX95AZGKdqzn/dvxxXsOd2OJnBZgQZ1dqtTpQSmQzBUkL8Re5OWFOKaDSfzl3fvrPq5RO2Kr+CMpuK2GdZXzXmKk34Lp4PrWvaRQgFAoKbKZ6rPocTm8uc7Ti/EsvM71AlVKJz0opUgKBRRLFPEM9yGVm5dmo0gKBRyf8DR+cAuM9JmlRCalFL914nkQEDzy+8fbfk1WcW+l/ZKJ5vWKhsJallI5lCi3M2iFZmbqrISQW9kPhJCbAfRm6nGDUMtsCdusNK7UiReGjT7Q56MZpHPFhiIpDJdl61bqhEIRFxbjNTN9dpO+64mH18pzHq/X2cBl80VEUjl87M17cPZ/vAU37OyXXYGQIw+pNZU6tbT1nvKFAQB3XD1Y93Eswx2s0Y7YKv5Iuqutl+1w8243XpyJrrt+rRiPKxHUGTadUEowLsBrX28lw5IeSpHOFVEsiVMy4dTWnSdXilNTIWg1BDftljeoG3KaEU7mkM0X8ezrS5gKJDG7nAalNSeeGrIyU998UDdoN0JDgCsx+VtB5UAJr8zNTjNB3YcAfJEQcpkQ8jqAvwfw68ouS12oZbZkMS5AryV1zccBSOqYG12SZyIpE4PNt1/GM3mUSu2fGNXKhSsJ5IsUB2r0nju6JK9dCRveZ21o1WCbxzG3BU6LHtscJmn2k9NbrLRfisbWG9nW2wqnfCFcvc3esC3JYdLBrNdWbUdsh5mltGQI3qscm/CgUKI4vUY1OJwUgy6lhFIS2cKmkuEPJLKrRFIAwGFWPulR+dobnYTdjDzuC+PQqEtqU5QL1lFzJZbFiWdmAIgBeieCJVLS3tn8MavTarDNYcL8cm8GdUomlzYrDYM6SunzlNKDAA4AOEQpPUQpfUH5pamHbvnRdEownsWg3dSwFYQFfRudTZ0OsKCuuUqdw6xHiQIJofcD7E544ORF/Oj8lVX3Tc5FAaBmULcRiYfpclBeL6hbXFM99tiMCCd4UNeLpARRFdKmokpdOlfAmcvLONZE+xQhRDavOqFQxEIs07PzdIzrx0UrmFO+0Kr7ldxMsZavjU4aykU2X0Q0na9jJaPc7xmveO1Ikgd1chJN5zA5F23q3NEqLKibnIviR+evSN6S8x2IpwSk8ZrmK3VsLZ28r5KwSh23NGiehkEdIcRLCPkygG9SSmOEkGsIIR/qwtpUQ5/FoIos2WI8u07uvhqs/XKjFcp8wQQ8NiP6GpilM1iFMbZJNgvVSOcK+OtHpvC3j06vuv/sXAxuq0HywFnLRrQIs0rrQiwDoVBdJn5ty8iA3YhUroi0gi1LnPZIlr041wd1vXu8PXNpCbliCcf3DjT1eK/DJEv75dxyBpQ2bwS+URh1Wty0242TU6uDOiXbnlY6QXr/mtkM7POy1iqoG0mPymNPDXsQNfHkdASUoulzRyuw6/T/96gP+SLF79+5FwCw0IF3XSAuwKzXSjPPzTLkMmOhR9svp4NJaDWkqlYApzrNtF/+M4BHAAyVf54C8HtKLUiNbHOYEIgLHfVDd4NAPNtUv3WvXHR9wWTTVTpgRcmplz1XOoVtUl+5Epey6YCY8Tsw4qyp0GXvspgPpRRTgQQcJh0oFTe51WCKrCzZwDJy4QTfoPQayTWVOkf5eItnejcAP+kLwajT4I07+pt6vNdhkqX9kgk29XpQB4gtmJcjacxUVNRDCQEagoat+u3QK0lDuZC6DZzVK3VK2h3FV7Vf8g4HOTk5FYLDpMOBDhQpa7HNaQIhosfgLXvcuLVcDZxfbl9YfjGeLb9ua6I8w31mLMay0mxmLzE5F8NVXjtMeu1GL0U1NBPUeSil3wJQAgBKaQGAPO6smwSvw4RMvtjzCpjBuNBkULfxQimUUkwHkk2LpAAr645mNm9AcLKiTeqJafH7lFDAdDCJ/XXk2rvdIhxMCEhkC3jTVaI4xUyNFsxgQoBRp5ECclYZ4HN1vUdtoZTe3ZyfnArhhp39TW8KWPtlpwk61nI81t/b7ZfASiXi1PTKuSWcFOBeI9EvF73S3i8XtQQqulOpW3ntCK/UyQalFKd8IdyyxwOdtpltcmsYdBqpsnvf0XG4rQYYdBosxNpPKInjNa1X1odcZuSLVKrO9wqUUkzORXFwlHvUtUIzn9YUIcSNshE5IeQogJiiq1IZzCIg2MOqfSmhgIRQaCqo02s1sJt0GzrzMLOURkIo8ErdGk75wjg24UGfRY9TU6Kq38sLcZQo6hp0iu2X3Us6+MrzkHe8QQzqallNLMbE6jHLLrIZnhCfq+s5UkIBhAAWgxggMRuDXk1mzUczuBhK4bYW2qe8DhOEQmlV9fGjD53B5x650NJ7+yNpWAxaVcyC7PJYMewySy2Y+WIJ5xdibW0Qm0FKvvVwUPePJy/hzr95HLkmTL2loM5eK6hT7nrEXtuk1/D2Sxm5HEljIZbFsQn5Wy8ZI30WeB1G3PkGLwghGO5wto1V6lpl2CU+ZyPn6r7x7Aze8oXHUagQT/JH0ohnC3W9RTnraSao+0MA3wewu+xP91UAv6voqlQGE3kIyCSFrQQr2cTmLtQbLTv97TNzIAR4cwMZ8kpcPW6k2SkL0Qymg0nctncAt04M4KQvLGWzANQ9+TlMeqRzxVUnTSXxlb2vbtrthsWghX+peqUuEM+usthglbpeyxpyRPVLq0EnBeBGnQZ6LelZoZQnylXtVjZm3goDckAUwXj01SCeu7zc0nvPLKUx7raqwgyaEILjez14ajqCQrGEz/zgFZyfj+PXb9mpyPutVOp68zydL5bwwKlLmA4m8cjLiw0fH4hnYdJrJLVLRjeSHuzYG++38qBORphy87VDDsXe49N3XYsH7jsiVQKHXCYstBlYUUoRaLITay1sXq3d95aDn74awFQgibNzK/Wisw3E3zjVaUr9EsBtAG4G8FEA11JKzyq9MDWxdiPQi7RqTNln0W/YRTdXKOEbz83ijqsHMdLX/EyKY5NX6phC3fG9Azg+4UE4KeDVKwmcnYthyGmqK2rAssbJLimDTgWScFn0GLAZMdZvqamAGUwIUqUbEJXxCOGVul4kJRRWGWkTQnrao/PkVBhehxF7W2jhZplulgS7sJhAoURb7sLwR1I9b2dQybGJASSEAv7k4XN48Gk/PnzrTrz3+hFF3sti0MKg02C5R4OQH78cQKjcFv7QaX/Dx7PN9NoAvhtJj0Q2D62GYLjPzNUvZWSmnIRUciZ237ATB0dXErFDTnPbgVUsk0euUGq7/RLAhtoaTJaDuUoV3nNzMRh1Guz1NmdpxRFpRv3yIoAPU0pfppSep5TmCSE/6MLaVIPXsXoj0Iu0Knfrshg2rD3mkZcXEU4KuOfoeEvPM+m1MOk1mzaoO+kTN6kTgzap+nDKF8K5uSj2N8hmdVt+fjqYwMSgDYQQjLstVdsvKaVS+yVDr9Wgz2KoW6mLJIWelWDezCSFguRRx+hVj85iieKJ6TCOTQy0VC1jLXQsQceq4IstzNkVSxSzSxlViKQwbtntgYYA3zozh9v2DuCTb3uDYu9FCCknDXszCHno9GWM9Jnx8TsnysbQibqPX6whQNaNpEciW4DdpIPbakCEC6XIhj+ShsOkk0TjusGQy4xgQmiq5XcttcR6msFh0sNu1G1YpW4xlkWwnMStVOGdnIvhmiEH9ArMNG5mmvlr5QG8mRDyFUII+4QPK7gm1WE2aOEw6VQR1DV70PdbN86m4aHTfoz2m3FbG/3sLvPGBaNKUixRPOFb2aRuc5pwldeOH0xeweVIumHfuaTE1oWqiqh8mcREOcM27rZidjmzzhQ+IRSQyRfXVY8HbMa6lbpP//sr+OhDZ+RfOKcuSaG4qlIH9K5H58VQErFMHjftcrf0vLXz0WdnxQxyNl9quo1uMZ5FrljCmIqCOqdFjxt29mP3gBV/+8HrFBFIqURs7++9z810MIHTl5Zwz43j+MAbx2DQaXCiQbUuWEdVWumkBwvq+m3i9brXFbjVgr/cPt1NhvvMoHRFEboVWCdWO+2X7L3nO7BT6ATWZnnrHg9emo0ilsmjWKI4vxDDQT5P1zLNBHVpSukvAXgVwClCyDjKoimcFbwOU08HdYvxLKwG7bpNWS1cFv2GzKZNBRJ49nXxotrIJL0aTrN+U1bqzs3HEMvkV3nmHJvw4Ny8uOlsdPJzdLFSF07mEMvkJZGbsX4LcoXSuvbkQIxVj1e3jHjs9St1s8tpvB5K8Q1Ml0kJBUn5kmE3dleAp1mYd9hIX2v+Ria9Fi6LXtokTc5FpQCn2fO7ZGegAuXLSr50/xvx779zqyQ4pSTi9aX3km8nTs/AoNXg/UdG0G814O37t+M7L8xLyq9roZSKlboabW9KJz0S2TzsRj3cVgPyRYpEl9rrNzszkVTXkzLMu66dLpRaYj3NMuRqv/WzU9g59qO37UKJAk9NhzEdTCKdK/J5ujZoJqgjAEAp/SsAfwLRs06ZZnsVs81p6mmhlGbtDBh9FgOSQqGtVoBOOHHaD4NOg/cfGW3r+c4NCkaV5tRUCISI2SxGZYC3v4GXDqvUdWMDzobMJwbFSt2OcsZz7VxdrTnPAZuxrqVBKCEglStuyuC9l6ka1PVo+2UoKW5y2jHP9tpFr7qkUMB0KIkbd4oed80Gdcy+Q03tl4DoP2gxtGZc3C69WKlL5wr4t+fn8Av7t8FdVuG99+g4kkIB331pvupz4tkCsvlSzQ4YpZMecVaps4rrXeJzdR1TKJYwt5zBji4fv50IltRKkDb/3qYNMyCfnIthr9eOo7vcsBl1OOkLNyX+xqlOM0Hd/2DfUEofBfDzAP5OsRWplEF7b1fqmjUeZ/RtgOdbSijgOy/M4x37t6Pf2l4v+2at1J30hbBvyLnq73LDzn4YdRrscFvgtNTPrnfTU8wXFO0MmMcg29yunaur5e/ksRkRTlRvJaJ0xU+Hz9V1l6RQgH1dUNebQinMvN7TTlDnNCEYz+Ll+RgoBd5yjRdA8+rG/qU09FqC7W3Mt2wV+qyGnhNK+d5LC0gIBdxXMct9eMyFa7Y78NDT/qrmzMEGs+pKJz3imTzsJrFSB3CvOjlYiGZRKNGuV9rZ+aKtoC6Rhcuib9uke8hlRjSdr1mRVgpRvTuGgyNO6LUa3LzbjZNTIUzOxWAz6rDLo65uh16gZlBHCLm6/O08IeQwuwFwA+BCKWvwOowIJoR1c0O9gjjM3fwGhw0Id7Pq9d2X5pEUCi0LpFTi2oRB3exSGmf8y7h9jb2DSa/Fr9w0jve/sXFVs5tCKb5gAnaTTlLi2u40Qach62wNFmsEdQN2IzL5IlK54rrXTgpiZhwQL76c7pFUVaVOgEGnWReENoPXbkQgLkiKbHdKQV1zn7dLoSRG+yyKmBZvFvosekQz+Z5poaaU4qGn/bh6mx3Xj/dJ9xMitoVdWEzgr6p4FUrnsBrJA4dZeaEUh0knJfu4rUHn+JfE5GO32y9NetHXsp2KWSAuNK1sXg3W+nmly9W6maU0Ypm8VJE7tncA89EM/vP8IvYNO9oawdnq1Lvi/QGAjwD4fJV/owBuV2RFKmWb04RiiSKcEjDYZl+zUlBKxfbLFjLHzEuoWxcJdlG9ZrsDh8faL7k7zZuv/fJrz8xAQwg+cMP64O1P335NU6+x0n7ZhUpdIIm9XrukOqjTajDSZ5ba0hjBeBYOkw5mw+rsIjMgDyeEdTOg4Yr2oo301dmKVGu/dJh0SOYKKJVoT12AwwkBAzZjWz5xXocJoaSAF2eXMewyY6TPAqdZ33RQ5wsmpSo1pzp9FgOKJYp4ttCVGb5GvDgbxStX4vjMu/et+8y869Awnru8hH94/BKu8tpx9+GV6ROphbxW+6XiQil5Uf3Sxq7XvTsCohb8G9g+PeRqT7AkEM82rWxe630BYG45gz2D3bMQYL50bHbu+IQ4XhJOCnjvYa7H2A41U4mU0o+Uv/0FSumbK28A3tad5akHFsgFe3CubjmdR65YammIts/KjLy7E9S9MLOMC4sJ3Ht0vCPDXpdFj0y+CKGwvsqjRoRCEd86I3r2bXe2JvpQiUGngVGnUdQIlzEdTEoiKYxxt1XKgDJqSYGzlrlqc3WVqpg8qOseQqGIfJHCtsbSwGHWg1Igmeutal0oKbTVegmI7ZfFEsWpqbA0q+p1GJsK6oRCEf5IWpon5VRnpROkNypLJ077YTVo8e7rqm8kP/XOa3HTLjc+8Z1zeHFmxYi+Vgs5w27SS0kPuaGUii3RJj3c5Zk63n7ZOTNLaRh0mrZFRzph2GXG/HJ1T9d6BOqI9TT7vkD3u18mZ6Mw6DS4atuKUjYLpvk8XXs00x/yVJP3bWlYa2MvztU1uvBUg1XqujXM/tDTftiNOrzr0FBHr+Msr3uztGD+57lFLKVyuO+m9ltSGd2Yf1pK5RBJ5bBnXVAnGpBXtlsF4kLVDPdARaVuLWyeTkP4TF03SQlikqRa+yUgzvb0EqFypa4d2OYoIRRwYJQFdSYsNpGwuxxOo1iivFLXgP5y0rAXgpDlVA4/mLyCuw+P1FSH1ms1+Pt7DsPrMOIjDz0vzR4F4lk4zbVnmRwmnWJJj1SuiBIVj0GzQQuzXsuFUmTgcjiFsX7LhnQeiCqUzXtiAqLdUShR/VraLIN2I7Qa0vVE6eRcDNdsX+1Fd6xcrePKl+1Rb6ZuGyHkegBmQsh1FXN1bwLQsC5NCBklhPyMEPIqIeRlQsjHy/f3E0J+Qgjxlb/2NXotNcAOqLWy7b3Aikdd85uclaBO+YtEJCngP84t4u7Dw+s2ja3CWnlim6QF88RpP3Z6rLhlt6fxgxvgMOkUr9S9HhZFUnYNrB5wHuu3IJEtrGqNDcazVVuVPXbxs1fN1oBV6vZ67bxS10XYJna9T133VFVbIZwUMGBvT2ypcnPErEK8DpMkilEPX3C18iunOkwR19fA2LsbfPv5WeQKJdzbYJa7z2rA37z/EEIJQVLDDDSYVVdylpkl6Bzla95GestuJmaW0l1XvmQMuczI5IstjZCEkwJKtLZYTzPotBpsc5i6ek1d8aJbHbx95NhufOqd17RsR8MRqVep+3kAfw3RvuBvIM7WfR7irN2fNPHaBQB/SCl9A4CjAH6bEHINgE8AeJRSOgHg0fLPqsdjM4KQ5hXSugkL6lqZ9TMbtDDqNF2ZT/vWmTnkio0vqs3gYkFdj1UO2uGVhTjO+Jdxz41jsmQNuyFqweYRxtYohzEj18tlBcxSiSKYEKpuiNxWIzQEVQ3Iw0kBGgJcO+Tklbouwj431czHK/+9FyiWKJZSufYrdRWbo30V7ZfNCGFNBZLQkPVJDc5qdritsBt10kzNRlEqUZw4PYMbdvRLLWD1ODLeJ6lhih519a2CWNJDiUo2O+bYMei2GRDmQV1HUEoxs5Red/3qFsMu8bPUyrVNStp3ENQBoq1BN6+pF0PMi251m+WY24Jfu2VnR2M4W5l6M3UPlufnfnXNTN1dlNLvNHphSukVSukL5e8TEM3LhwG8C8CD5Yc9CODdHf8WPYBeq4Hbamwqm1ssUUTq+HDJDQs0W/Uw6bMon/krlij+5Vk/btzZjwlv59ltVqnbDGIpJ57xw6jT4H3Xy2ML2Y32S38kDUKA0f7VWTbWJz9TVsCMpHIolGjVlhGthqDfakCoSitROCmg32rEaL8ZwYTQdR/FrUqq3D62vv2yewI8zRJJiZnrdmfq3FaDGJh5rNL5ZJtjRQirHtPBBMb6LW1Li28VNBqC/SNOnNvgoO6kL4SZpTTuOTrW1OMJIbj36DguLCbwwswygg2sgrpRqWPHoFip672kspoIJQWkc8UN85hsx6tuMcbGa9qfqWPv3U2vurOzohfdwVHeZikn9dov7y1/u4MQ8gdrb628CSFkB4DrADwDwEspvQKIgR+AwdrPVBfbnMam2i+//MQlvOlzjyGb746YRyCeRb/VAKOutY1Gn9Wg+CD7izPLmF3K4JdvbO6i2giXZXNU6iil+PeXFvD2A9slUYFO6UalbmYpjSGned1nbazfAg0Bzs+Lm7hG1WOPzVi1UhdKCBiwGzHkMoPS3pxh3YwkhVpBXe9V6phHXbuVOp1Wg9F+yyppe9ba1EgIyxdIdlU9Ts0cGHHhwmJ8Q0WtvvviPPqtBrx137amn/OuQ0OwG3V48Cl/zW4DhpJJj/iaSl2/1cBn6jqEKTR3286A0U5QdyUmV6XOjMVYtqofoxKcvrQEu1GHnR4+fywn9dovWf3ZBsBe5dYUhBAbgH8D8HuU0ngLz/sIIeQMIeRMKBRq9mkbitduaqr98scvB5AQCphdal3lqB0C8azkGdYKfRa94kIprF2Pqcx1ilSpU3lQt5TKISEUZPu7ACyoU/r/UxwyX4tJr8Udb/DiOy/MQygUK+Y8q1+IBuzG6jN1yRw8NoOk1sVbMLtD7Zm67pnaNwv73LRbqQOAEx+6Ef/vO1bsQlg1hmXFq5EvlvB6OIW9XCSlKQ6OOJEvUrx6ZePm6i4sJnBo1NVSwtNq1OHuw8P4weQCiiVadzOtbKVOfE0Ha7+0GhBJ5XrG+0+NXGZ2BlWuYd3AbTXAqNO0dF27GErCZtRhoIPzHSAqYOaLtGoyVW5i6Tx+eG4B7zg4BG0PWeFsBuq1X/5D+eunq92aeXFCiB5iQPe1ipbNACFke/nftwMI1nj/ByilRyilRwYGBlr5nTYMr9PUsHIQz+bxYrns7I90K6hrTxmpz2JQXCiFZaRYhqpT7CY9CAFiPSKV3S5MWliuvwsAOEz6rszU1Wpdue/oOCKpHH50frHCeLz6hWigRqUuXFGpA7itQbeQgjrTWp+68sxQD1Xq2Oem3UodAIz2W1b5p7GNeyBR+/zuj6RQ4MqXTXNgVJylmZyLbsj7F4olXAql1tmvNMO9R8fBChr1BCqUTHqsbb9024wQCiWkc5vDzmcjmImkoCHASN/GBHWEEAyXFTCbRewOsHU8g9bNROm/vjCHbL6Ee5tse+Y0T0NLA0LI31a5/S9CyLsaPI8A+DKAVymlf1PxT98HcH/5+/sBfK/dxfcaXrsJS6lc3XaSp6YjUnnb36VK3WI825bnisuivJH3QiwDj80g2wyKVkPgMOlV337JTqzDMgZ1dpMe6VwRhaIyc2hJoYBIKlezdeXWPR6Muy04cdqPQFwAIStG42vxlCt1lVlnSilCSVGqfns5ScGDuu4gCaUYVgd1Rp0Gei3prfZLGSp1a/HYxDm7ep0YUwFR+ZUrXzbHkNMEt9WAyQ7n6h6fCrXV9TKzlEauWFpnv9IME147ju7qB1DfKkjJpEc8s779EgBXwOwA/1Ia251mGHTNuH0pg2hAvv66thDN4KevBNbd7wsmZekOGOpSUEcpxddO+3F4zIVrh/g8ndw088k1ATgEwFe+HQDQD+BDhJD/Xed5twC4D8DthJCXyre3AfgsgJ8jhPgA/Fz5500BqzrUK1+f8oVgNWhhM+owE0nVfJxcXA6nEEoI2NmGGlufRZypU8I4lTEfzcpajQLEFky1t18uKBLUiRd/Nh8lN/7y53m8hnKYRkNw743jeO7yMk5OheCxGVf501TisRkgFEpIVKw1ni0gVyhhwG6ESa+Fx2bg7ZddYsWnbnXyhRDSFQGeVgglBJj0GlgN8omV6LQaeGxGBOq0X/oCSRAC7B7glbpmIITgwIizo0pdrlDCb3z1DL74s+mWn+sLloPwNgW6futNezDkNGGnp/a1VcmkRyKbh1ZDYC4nRN3loK4XvP/Uij+Sxg7PxlTpGDs8FkwHk8ivSb5+4SdT+I2Hzqw61y6ncggnBVkSSWP9FthNOvzD4xeRVsBXkfHUxQguhVOyqJ1z1tNMULcHwO2U0v9DKf0/AO4E8AYA7wHwllpPopQ+QSkllNIDlNJD5dt/UEojlNI7KKUT5a9L8vwqG4+3XD2o1YJJKcVJXwg37XaLZsxdqNR97Rk/dBqC91w33PJz+6wGlKjYMqoUC9EMhpzyBnUui/ordQvRDMx6rST8IgdKi1qwIfN6ymHvu34EBp0GL81G6woMsPmASgNyqQJTru4Nu8yYb6FNYc07ywAAIABJREFUhdM+qVwBRp0GuipBeDcEeFpB9Kgzyi6J7XWY6rZf+oIJjPZZYJYxmNzsHBhxYTqYlNp7W+XCYhy5QkmySmmF6XJQ106lDgBu2zuApz55x6o23bUomfRIZAuwm3TS53ylUscVMNtlI+0MGLfu8SApFPDS7Eqyg+0dKQXOz69IU7DExB4ZKnVmgxZ/+4Hr8MqVOP7425OKzWY+9LQffRY93rZ/uyKvv9VpJqgbxopoCsrfD1FKiwD42aMC1uJYq0XHH0ljdimD43sHxKBO4Zm6bL6Ib52Zw1uu9dZtEalFXzmgUEoshVIqBnVKVOpUbmmwEMtgyGWSdWMqeSYpFKSzJEW9oK7PasA7DwwBqK/WxQK3yqq3NCtVDviGXGbeftklkkJhnUgKoxsCPK0QSgo123o7wesw1m2/nA4m25rP2socHHWiRFdUcVuF+dzNtHEtnQokMOwy1/xcy4VDoaRHIpuXEnWA6O8JABGugNkW8WweS6nchtkZMG7a7YGGAKemVgQCfcGkdO6prGz7gqLIkFznnTdfPYhP/sLV+OG5K/jbR1uvfjdiMZbFT14N4P1HRrnti0I0E9T9FYCXCCFfIYT8M4AXAfw1IcQK4KdKLk5tsMpDrUrdKZ94kB6bGMBYvxVzy2lF5WN/MHkFsUy+7TJ3X1lKXymxlGg6j3SuiCFXZ1K8a3Ga1V+pm1+WP9hlKmlsFkNu/JE0+q0GKXisxX03iZ/HegIDUqWuYoOytlLHgjqu9qY8KaGwTiSF0Q0BnlYIJ9o3Hq+H11FbCIuJbsiRMd9KMOPhc20GdZPlasaVeLZliyAmMKE0ilbqjCvn2n4bb7/shJkNVr5kOM16HBp14XFfWLrvZDnAc5h0q2ZQfYEkrAatrGMav3FsF+4+PIwv/HQKP3utqo5h23zlqddRolQ2CyvOehoGdZTSLwO4GcB3y7dbKaVfopSmKKV/rPQC1US/1QC9ltT0qnt8KozRfjN2uC0Yd1uQL1JFKw0PnfZj94AVN+1yt/V81vqnlFcdm4ca6ZO/Uqf6oC6alfVEDShvFD2zVN3OYC0HR5z42Jv34O46LcEscKu0NWCVOk958zLkMiOdK6r+/1oNJLMFWA31KnW9E9SFkoKsIikMr6O2EBYT3eAiKa3hsRkx7DJLFbdWmZyLQashoBSYW26+WlcsUVwMdaeyqtTxwdovGVaDFgadhgultMmZy+IkUC+IdxzfO4DJuai09zrpC2P3gBXHJgYwOb9SqZsOyqN8WQkhBH95936Muy34+zZmVWvx01cCeODkJbznumGMuze2xXUzU898/Ory18MAtgOYBTADYFv5Ps4aCCEYtJuqGtTmiyU8fTGMYxMDIIRIJf4Zhebqzs3FcHY2inuPjrd9wLNK3VJKmU2z3HYGDFG1U1mBFyXJ5osIJwXZ/y5Kz9TVszOohBCCP/r5q3BkR3/Nx/RZDNBqyKr2y3BSgFZDpM/lcLnCy8VSlKd++6Ve0bnbVsgXS1hOK1WpE1+z2vmdKV9yj7rW2T/cnlhKOleAL5jAzbvFpGUr4wxzy2kIhRL2timS0gpKBXXxbH5VVwQhRPSq4+2XbXHKF8a427JhxuOVHJsYAKXAk9MRZPNFPHMpgmMTA9g/4sTsUkYK3KcCCexRIJFk1GklUbMLi03bS9fktcUEPv6NF7FvyIm/ePd+GVbIqUW9St0flr9+vsrtrxVel2rxOoxVDWpfnIkilSvi+ITouccyFUrN1Z047YdZr8Xdh0fafg22eVaqUqdYUGcWBV6SCio4KQn7/MhfqWveM6mZgLiydThXKGEhmpGtdUWrIei3GlZV6sKJHNxWAzRls1JJgnmZB3VKk8oV1ilfMnqpUreUyoFSee0MGGwuOVhFLGW6PNvClS9b58CoE/5IGv5ICkupXNOdBC8vxFGiwDsPijO6rVxLfQH5BCYaoWT7pWNNS7TbZuBCKW2QK5Tw9KWItD/baA6OOGE36XDKF8KZy8sQCiXctncAB0bEKuLkXBSxdB7BhKCYL+b7rh+BUafBidP+jl5nKZXDh7/6HCxGHf7xV45wISmFqWc+/hvlr2+ucru9e0tUF9udZvgjqXVzPo9eCECrIbipnFXc5jDBoNXAvyS/rYFQKOJ7Z+dx18GhuspcjbCbdNBqiGIzdQuxLAw6jSTFLBfsd46pVCxFqWB3pf2y/gb8R+cXcejPf1xzfkh8zBVc9+c/xvP+ZQBi5rtEIWtbxaDdiIWKBEmorGrI4Abk3SMlFGGtU6lLCgVF54ObRQ7j8VqwoK6aWMq5+RiGXeaafyNObQ6VTchv+9xjOPy/foL9//PHeHI63OBZwNnyPN2brhoQLYJa6Hrxdah82QrKtV+uFkoBgH6rcdUcciO+++I8bv7LR9fJ5281nvcvI50r4tiEZ6OXAkC0ULl1jwcnp0I46QvBoNXgxl392D/sBCFi2/F0SEwkKdUd0Gc14B0HhvDwC/Md2SB95gevIBAX8MB912ObU179BM566rVf3l3v1s1Fqok73jCIhVgWT12MSPcJhSL+9cwcbr96UAo4tBqCkX4z/GH5K3WXQilk8yXc0uEJSqMhcJn1iqlfzkczGHaZZZced5ZnAdU6a6WE8TgAGHQaGHWaVd5v1TjpCyGeLeDrz85U/feXF2L4/W+eRTxbwD898TqA5pQvW+XQqAvPX16SNhzhNaqGbqsBRp1mVeDHUYakUFi3gWQ4FPY/bIVQkimkypsoAlaCurWdGKGEgP+6EMRb922T/T23Akd3uvH5XzyIT991Lf7HO64BAJxtoh1zci6G7U4TBu0mjPVbWrI18AUS2OYwSebgSmI36ZHMFWQdB6CUIikU4FiTtB12mVuaLfQFE1iIZVd1RGxFTvlC0FUk3XuBYxMDWIhl8a/Pz+HIjj5YDDrYTXrs8lgxOReTqs1KzvHed9M4UrkiHn5xvq3nF0sUP301gHcdHMJ1Y30yr45TjXrtl+8s3z4E4MsA7infvgTgXuWXpk7etn87+iz6VSXrH51fRCSVW6dCOd6vjFedZKoqQxaSzacpgWhnIH/mRqrUqTSoW4hmQQjgdcpfbWimFehcWbTg68/OrMvghpMCPvLV5+Gy6HH34WE88vIiAvGspBwm5zzCsYkBpHJFvFCuBoYSqyt1hJCyVx2v1ClNI6EUQDkBnlYIS5U6+c8rfRY9DFrNOq+6b52ZRb5IcQ9XdGsLjYbgvdeP4P6bd+DXb90Jj83QlEXBufmY1I427ra0ZGvgCyYVa1tbi8OkA6VomExrhVSuiBLFukTLuNuC5XS+6RnXdE4U/ak2MrKVOOkL4fBYX0Pl5m7CqoZLqRyOVbSFHhhxYXIuCl8wCZNeI3vyt5KDI07sG3bgxNP+tlSmJ+eiiGcLOL63N9patwL12i9/jVL6awAogGsope+llL4XwLVdW50KMem1+MUjo/jxKwHpRHnitB/jbguO7VldORt3WzFTpVWzU6YDCWgIsNPTeStcn8WgmJrW/HJGkRPSimqneGGbCiQw2wWjd7mYj6YxYDPCqJO/99xh0iFepxVIKBRxYTGOq7fZEYgL+OkrgVX/9psPPY9ISsA//soR/O7tEyiUKL7x7Cz8kTQsBq2sbW8373FDqyE45QuDUrquUgeo16vucjjVUkZ9IymWKDL5+u2XgHICPK3AWs88ClTqCCEYdBhXCaUUSxT/8swMbtnjxi4+TycLY/2NPVxjmTxeD6ckS4QxtwWzTVoElUq07CnYHaVSJZIe7LXWBiE7mABbkwFuWhCDunr+i5udSFLA+fl4z7ReMkb7LdhV3sNVru3AiBPBhIAnfGHsGbRJM+ZKQAjBfUfH8VoggTPl5GornJwKgxDRUJ3THZrxqdtBKb1S8XMAwF6F1rMpuOfGMRRLFN94bgYXFuN47vIy7rlxbN3BN9ZvQSpXlN1XxhdMYtxtlcXc0eswYXZJfi8woVBEMCG/wiOwUqmLZnK4sBjHu7/4JD71/Zdlfx+lWIhmFfm7AI3nO169kkC+SPGx2/dg2GXGiWfEijOlFH/23fM441/G53/xEPYNO7HDY8WxCQ++/uwMLoWTGOu3yNpK6zDpcd2oCyd9IcQyeeSLdFWlDhDbjdQUsDN+5+sv4hP/dm6jl9EUqbLgUD3zcaA3grpQQoDVoIWlRlWxU3YN2PDYa0HpM/ezC0HMRzO4r00vUM56xt3WhvNxrJuAVep2uK3IFymuxBoneOajGWTyxa5V6pRIejCv0bWVurH+1gTY2LFdb356s/NEeX6zF6tJb923DeNuC67Z7pDuY4mM1wIJ7O1CYuKug8Mw67X4j3NXGj94Dad8IRwYdqJPZt0ETm2aCeoeI4Q8Qgj5VULI/QB+COBnCq9L1Yy7rTi+dwBff3YGDz51GQadBr94/ei6x+3wiFk1uRUwfUH5/Hdu3NWP+WgGl2VeYyAmZgaVCF5cZvEE8noohQ8/eAbpXFEx6wglWIgqU8EEGrdfMmnx68b68Ms3juHJ6QguhpL4pycv41tn5vC7d0zg7Qe2S4+/7+g4FuNZnJwKyTpPxzg2MYBz8zFJMp551DH2DNoQTuawrCJvpkKxhNcCCVmkortBqtw21rhS1wPtlwp51DE+fde1KFHgww+eQVIo4KHTfngdRtz5Bq9i77nVGOu3YCGWqeoHyGAzdweGxQ0uU91tpkI1LeN4QjMokfSoValj7e/Nzhey9sutHNSdnArDZdFj3/DG+9Ot5Q/fchUe+b3jqwoC1w45oC3/3A31VrNBi71emzTD1yzxbB4vzkZXtY5ylKcZ8/GPAfi/AA4COATgAUrp7yi9MLVz39FxBOICvv7sLN55YKhqpoJl1WZkVMDMFUq4HE7JloVkEr+nfCFZXo+hlBgIAJj0Ghh0Gnz5ydcRSgg4uqsfC1H5q41KQCkVBWRkNmRnNKrUnZ2NwWMzYMhpwvuPjEKvJfjkv53DX/zwFfzCvm34vTsmVj3+9qsHMeQ0ya58yTi21wNKge++JA5qr63UsYsamyNVA7PLGeQKJYSTOVUYBbOgzlZDKKXXKnVrW3TlZKfHii/+8mFMh5L49X9+Did9IXzwhjHotM3kRznNsMNjAaXA7FLtqtvkXBQ73BZJFGslmGkc1PnK9hPda7+UP+nBjrW1lTqbUdf0TCIgev0BW7f9klKKU74Qbt3jkQKlXkKrIes6rkx6reSv2K3P8J5Bu3TcNMvTFyMolmjPtbVudpq6ElFKH6aU/n759rDSi9oM/P/t3Xd82/Wd+PHXR8OWt+MZ78SJkxCyYzLIAFoopaWMMhMIaa+F0tLd3l3HtRd6d7R3ba8tLS1lFQgQflBWj5YWykoocXbibDvLjp3EM962bEmf3x8a8ZBsyZYsy34/Hw8esSVZ+sjI+ur9/byH+8MuwJ3LvBfQ56XEoFRwd+pONbRjc+ig/bEXpMaSlxLD5rKhW0wHIlRt+8GZB54UY0Zr+J+b53HlRZl0dNsjonFKY3s3VpvD89oJtkSLs/GNrwC3tKqJebnJKKVIT4jm43Oy2H6qkZmTE/n5rfMHpBCbjAbWLHG+vvODNKOut3k5SSRaTLy+7wwwsFW9++AW6AEnnMprLqz1WAQEo+4PkPGDzKlz3i78f1/1bdaQjDPobWVRGj+8djbbTzZiUIrbL5EGKcHkz8nO/VXNzHWloYFzlJDZqPwaEVRW00ZGQrQnIAy1xBCc9HA3Quk/pw5cNYl+niieKDt1p+rb+ckbR7D1a/x1tKaV2lbrmJlP56/5rrTj0dptLsqMp6bFGtBnqM1ldcRFGaXr5SiT04shYjQovn31TG4tzvXM4ekv2mQkK9ES1KDOM1Q1SH/sSilWF6Wz9Xh9UGfZuIO6rBAFL9fPz+a718zi+gU5nt3AqggYUn2myXlwDVVN3eKCSdS3dfPg28cGXNdutXGsro25vdJQvnzFdD46K4NH71rss05p7dJ8Vs9ID0kxtMloYGVRmqe5y4BGKUkW4qKMAaeGhFPvXcWymrEfjLa7min46n7pbgs/WAOe0VLXZg1Jk5T+7lpewDeunMFXPjJdZi8FmTuN29dxsa7VypnmLs8HW3Aeb/Mm+dcB89CZFmZOHp0dDgj1Tt3AwHRKapzfO3XuXfjxHtQ9vbWCh98/zt8P1/a5/JXd1RgNistnRlZQd+PCHD45L4u8EJxI9cY9C+9YACdPt5TXs3xaGlEmCTNGk/y2Q+jTi3L5n5vnD9o8Ij81looA5usMpby2FaVgWhA7sfVvLR8MZ5o7SYuPDkozF2/+7drZfOGyaQCeVMZI6JJYHcIdTIBbinP59KIcfvH3Mt7oV/h8oLoZrWF+3oUPSzMnJ/D4Zy4hd5Lvg0dqfDRP/9MSpgSh26o37px8s1F5muC4KaWYnhEfETtebuU1rWS5gtFIWHfbEDV10SYDUUZD2NMvu20Omjp6QjLOoD+lFF+7soivXyk9w4ItNS6KuCijz6DOXfc7L7fvyVLnsXTwYKarx87Rmlbm53o/0RoK7p3sYJ708JV+Cc7fw9mWrkFrEt0myk6du3zk2W0XRk119dh5Yedprrook4zEyDoxs7QwlYfWLhq1lFF35pe/J08rGtqpbOxg9QxJvRxtEtSFWUHK0J2+AlFe6+xCGBMVvGCpd2v5YKlu6iInBDPqvHEHSJEU1IWqUYpSigdunMvC/GS++cI+Dp5p9lxX6uooNzdn9D7w+MOdk58aF+21ffNw8v3DyTkjK4HpmZGxbk9NnY+gTinlqtUMb/plQ7uzLmg0dupE6CilyB+kA+a+qmYMytkworcprp8ZrHb60NkW7A7N3NzRa4phMRuDftKjtasHo0ER4+WkaEHq0DWJbu6grqXLRmf30EFgJDrT1El5bRs5yTFsKa/nZL3zJPpf9p/lfEcP65ZL59qh5CTHYDEb/K5d31zmDKKlScroCzioU0o9pZT6nVJqTigWNNEUpMVS39YdtA9Ex2qC1/nSzd1aPpjNUqrPd4RsN6q/1LgookwGzkTAgNUzTZ3EmI2eWXuhYDEb+f26xSTHmrnn6V3Utzk/DO+raiI7yTKgGUm45U6KpTA9joxE7+uaMYx8/3Cxe2ZkxVOUEXhHsXDwjDTw0SgFnDsG4U6/dKcuh7qmToReQUqszw6O+6uamJ4RP2DnOD8lljarbdDmQ6Wnnbt8o7lTBwT9pEdrl40Ei8lrFlAgDdg6um1kut5Xx+tu3Qeuk9H/fdM8TAbFsyXO3bpnSiooTIvj0mmp4VxeRDAYnBkx/pQL2B2a1/aeIXdSjGduohg9w9mp+w3wd2BdkNcyIc1y5fYfOjPy9uY9dgcn6tuYHoKOSKuK0imtbg5Ktz6tdUhnsfWnlCInOcazCzaWnWnqJDvZEtR5b95kJFh49K5iGtqt3LtxF1abndKq5gEpTWPFf980j+9/4iKv1xUNI98/XKrOd2C1OZiRGc+MzHhqW600d4ztYLRtiJ06cKbg1reGt4Pea3uriTIaKJ6SEtZ1iJErSIulqrFzwDBxrbXP96kCPzpgllY1k5EQPep1kEN1HQ5Ua1eP19RL6PV7qB88A6jb5qDHrilMc75/nhunQd3m8joyE6NZMT2Vq+dM5sVdVeyqOM/uyibuWFYQ8mPteFGUkeBXucBP3jjMzorzfPmK6fK7DYOAgzqt9Q6t9Uta638NxYImGneqmzv1bSQqGjroseuQdERa7Wot/49jI0/BbOroobPHPmpBHUB2siUi0i/PNHWSM0j9WjDNyUniZ7fMZ2fFeb75//ZR2djBvLyxN6sH4JIpKSwt9H5GNdB8/3C60Mgo4cK6x3gw2tZlw2hQRA9S8D450UJNa/g+FLZbbby8u5pPzssiRQbdRryClDi67Y4BgUZ1UycN7d19mqR4fsYVzAy2Q7XP1d13tA01HzRQrV02EqK9Z3O4axKHKutwp1u6a6HH406d3aH54Fg9q4rSUUqxblkBzZ09fOW53VjMBm5elBvuJUaM6RnxnG3uGvR1/MddVTy65STrlxdw+xLpChwOvk+9DkIp9YjW+p5gL2YiSk+IJjvJQmn1yIM6906Fu817MM3LTSYpxsyW8jo+NT+7z3VWm529lU04/BwDd9p1sAlV3Zg32UkxvF8W3Fl7oVDd1MXsfrUioXTtvGzKzrXy4DvObpijnZYUDIHm+4eTe43TM+JpcaWLlte2jendpXarjbgo46BnXTMSo3nv6Oh8KHQ4NKca2ins1Qzq1b3VtFlt3LlM6mPGgwsdMNv7HCfcJz+9BWa5k2IHHRHU2tXDifp2rl+QE4IVDy74O3U2nzt1SikKUuOGbMDmTqsudAV1teNwVt2B6maaOno8ddlLp6Y4095r27i1OHfUxlqMB+7Plcdq27yOKdhV0cj3Xt7Piump/ODa2aO9POHiM6hTSvn6lKGAT4RmORPTvNxkT0evkXDvAkzLCH4XQqNBsXpGOm/sP8cXLpvm6a7ZY3fw2T/s4MPjDQHfZ2F6aLolepMzKYbaVitWm51oU2g6bo5Uu9VGfZt1VINdgK9fOYOjNa28e7SOOTljc6duMIHk+4dbeW0rmYnRJMWYSYg2EWMe++MY2qz2QVMvwblT195tp81qG/K2I/XOkVo+//ROfvzpuaxZko/Wmo1bK7goK5FF+ZF3UkIM5J55WdnQwaXTLlxeWtWM2aiYlTXwxKXFbCQ7Kcbn39OB6ha0hnmj2CTFLcFiGjId0l92h6aysYNFBb5f6wWpsRwd4v3QPXg8IzGa2CjjuEy/3FJeh1J4xu0opVh/6RR+8NoB1i2bEt7FRRh3Blh5zcCg7kxTJ1/YuJvsZAsPrV2EySg9GMNlsKNvHVCBM4hz067vM0K5qIlmXl4Sfz14jqaObpJjh586VF7bRu6kGJ/zxEbqX66eyT+O1XP3Uzt55UsrSIo18x+vH+LD4w187xOzAuqamGAxhWRH0Rd3qmdNs5X8MVq8e8C1W3tx9uh+6DAYFL9Zu4gzTZ0DRgZEihkZCWw9EfiJhdF2rLbN87p3B6NjPf2y3WobtEkKQKarJfi55q6gzcj0xd1A4wevHqAwLQ6TUXHkXCsP3DhXajjGiexk5zDx/vVxpVVNzJqc6PPE3NLCFN49UovdoQe0e/c1CmE0JFjMnoHhI7W5rI5zLV18Ym6Wz9vkp8by9mHvvwe33vMnMxMt4zL9cnNZPXOyk0jt1Txp7ZJ8VkxPY2qIRvCMV3kpsUSZDAOOVx3dNu5+eifWHjvP37N0RJ9hxcgNdqQ+AXxUa13Z/wql1OnQLWnicae8lVY1s3rG8FvAltcGv/Nlb3kpsTx852LueKyEL2/azVWzM3l6awV3r5rKPaunDX0HYeTe/apu6hyzQZ1npEAYziSbjQYKUiP3IDc9M56X91S7GgiMzcDU4dCU17Rx+5I8z2VFmfF8eGxsB6Pt3TafM+rc3J1Ja1tCH9TVtlqJMhnInRTDF5/dzcXZicRHm7h+QfbQPywigtGgyJ0U26c+zuHQ7K9q5rpB/j+vLkrn5d3VHDwzsJlKaVUzeSkxYam5DGb65caSCtLio/nY7Mk+b9O7JtFX5od7nEFstJGMhOhxl37Z2tXD7srz3LO6sM/lBoOSgG4YjAbFtPT4PmUODofm2y/u49DZFp5Yf0lImvSJwAy2R/pLYGDirNP/hGAtE5Y75W3/COrq7A7N8Trn/KtQWjI1hf+4fg5byuv54WsHuXxmOt+5xntXwrEkHLPqunrsbCyp8PsM7b6qJnKSY0iTluwBczcdGcvDvKubOunssXvWCs51n2vpCtpZ/FBo7Ro6pXKya6duNJqlnGvuYnKihcfuKsZmd7ClvJ6bFuUMGXiKyJKf0neY+MmGdlqttkHrfle6aqe8zVTdV9XEvDDN4EywmGmz2gZ08wzU6cYO3j1ay5oleUQN0rjIU5NY77uuzp1+GRdlYnKSZdylX2493oDNoWVWWhD1H8Pz4Dvl/GX/Ob57zSyumCUJfGOBz3cFrfVDWut9Pq77deiWNPEkxZiZmhbHvtPDr6v78Hg93TbHqNQL3L4kny9fMZ1LpkziwTULfaZ3jCVZrhbWoxnUPb+9kh+8eoCvPLfHr4P5/urmsNR7jAeefP8xHNS5A073CAa4sO6xHIw6G6X4m34Z+rP9NS3OoK4wPZ7f3rGYi7IS+cyKqSF/XDG6pqTGUtlwYZi4O31ysEyGtPhoLs5OHNAUq7G9m6rznWF7f010pS+7x4MM13PbK1HAmiE6C7prEisG6YDZ7tqpi4s2etIvBxvcHmm2n2wk2mRgcYGvvQkRqKKMeKqbOmm32nhj/1l++fdyblqUy92rCof+YTEqfAZ1Sql/6fX1Lf2ueyCUi5qI5uUmjWiswTMlFaTERXHV7Mwgrsq3b189kxfvvZTEMZrq1p/FbCQtPmrUZtVprXlmWyXJsWbeL6vjJ28cHvT2TR3dVDR0jNk5cWNdXkos0SYD5WO4WYq7FqF3irRnxt4YbpbSbh06/TIu2kRCtGlU6nJqWro86Z4ri9J442urJJ1qHMpPjaPVauO8a45jaVUzFrNhyBKD1TPS2V1xvk8AFc56OsBznBzJWAOrzc4LO07z0YsyhxwH5K5J9NUJFKDD9fuJiTKRkRCN1eaguXPsZgwEqrSqmdnZiYPuaIrAuDPB/rTvDN98YR+L8pN54NNzpJZ5DBns1X57r6+/2++6j4dgLRPavNxkzrV0UTuMD0Vnmzt561ANtxbnjdnOjmPBaA4gLznRyLHaNv7tk7NZv7yAR7ec5MWdvktR3QG9t/lLYmje8v3HmrKaNtITovsUkudOcgajY7lzZ5vVd/v03jISo6kNcfql1pqaFqsn3VOMXwUp7mHizhTC0qpm5mQnDdlZb1VRGjaHZmuvjsylVc0oBXNyRm9cTG/uv5+R1NX99cA5Gtq7WefH2A6jQZHAQ35rAAAgAElEQVTXryaxP3dNXVyU0TOMvWac1NXZHZoDZ5ojckTPWOY+Cfn9V/aTHGvm4XWL5TPnGDPYu6Py8bW378UIuT/MD2e3btO2SjRwx1IZ9jiY7OSYUUu/fKakgqQYM9fOy+IH185mxfRUvv/KAZ8f3t1nki+OwJECY0VRZjxHz7WOyRQiu0NTWtU0YJfB6OmAOTaDUa017d124qKHPnBnJlo41xzaoK7VaqOzx+5J9xTjl7su7EvP7OajP3+Pvaeb/GoitbhgEjFmI1vKnSmYWms+PF5PYVpc2JooJXh26oYf1D1bUklBaqynPf9Q8lNjB9+pc9XUxbq6X8L4GUB+vK6Njm67lDMEWUFKLGajIspk4NG7islIkPfhsWawoE77+Nrb9wMopZ5QStUqpQ70uixFKfWWUqrc9a8kO7tcnJ2EQRHwvLoeu4NNO05z+Yx08lLGZlfHscIZ1IW+bqC2pYu/HTzHrcW5WMxGTEYDv7p9ITaHg9dLz3r9mdKqZgrT4iJ2pMBYsGJaGmebu9hx6ny4lzLAj/9ymLKaNm5alDvguqKM+DFbU1fXZsXu0KT70bxncqIl5Gf6a1xBozv9UoxfhenxfHbFFBZPmcSsrEQ+OTeL2y8Z+sRltMnI8mmpnmYpv33vOCUnGv362VC5sFM3vPTGpo5udlQ0csOCHAx+1rBnJliob/P999jebcdkcH5Az3R9OB8vzVLc/QmknCG4TEYD37nmIh5ZVxyRM20ngsGCuvlKqRalVCswz/W1+/u5ftz3kwxM0/wO8LbWugh42/W9AGKijMzITGBfgDt1bx6soa7VyrrlQ6dkTHTZyTF09thp6ght3cDzO05jc2jWLr3w/yQtPpp5ucmes8f9lVZJk5SR+tT8bBIsJp4pqQj3Uvp4cedpHvvgJJ+5dAo3LfYS1GUmUN3UOeImCqFQ6TrTX+BHzVpGooXa1i4cI+zwNxh30Cjpl+Of0aD4909dzENrF/HQ2kU8uGYhMyf71915VVEaJ+vbeeKDk/zszaNcvyCbz68KXzOdkaZffnCsHq0JaORRcqx50GNdZ7ed2CjnDnzvkSTjQWlVM/HRJgql1jboPrdy6ohGb4nQGqz7pVFrnai1TtBam1xfu78fcjtBa70ZaOx38fXAU66vnwJuGPbKx6F5uUnsr272upPU1WP3+jPPlFSQOymGy2ZIO9mh5CQ7Pwh6q6uz2R1ef8cOh6az2/vv3hub3cGm7ZWsKho43HT1jHT2nW6iud+Btrali3MtXXJWcYRioozcvDiXNw6cpa51bNSG7Kpo5PuvHGDl9DT+7ZPeR3+457odH4O7de70rQI/sgAmJ0bTY9ec7+gO2Xrc6WGSfikG4/7Q+aPXDzEvJ4n/vmleWJs5JIywUcqWsnoSLaaAaq4TY8xYbd6Pa9C3AZLFbCQ51jxuaupKq5qYk5Po966mEOPFaLcFytRanwVw/SuRSC8L8ibR2N7Ne/3aMb+w8zTz73+THaf6xsjHalvZeqKBtUvzI2KsQLhl9xpA3t/3XtnPDQ/9Y8DlD28+zuqfvku3zeHXY2w5Vs/Z5i7u9FLMvrooDYeGfxzvO0PJvTsrO3Ujd+eyAnrsmhcGaUozmv75xVKyky38Zu1Cnw0eZrg6io3FZikVDe0YlLOhy1Au1OWE7oPhOQnqhB8K0+LInRRDRkI0v19XjMUc3mYO7p26lmHs1Gmt2VJex4rpaUM2iektOdYZSPrqaNnRa6cOnLvf4yH9stvm4PDZVmmSIiakMdvrVSl1j1Jqp1JqZ12d95S18eaGhdnMmpzAVzft4Xid86z9jlONfP+V/VhtDp744GSf2z9TUkmU0cCtxXnhWG7EyfExgLy+zcore6o5cq6Vioa+3cL+euAcda1WdlX4V6f1/tE6LGYDl3lJT5ifl0xCtGlACmZpVRNGg+LibAnqRmpaejyXTkvluW2VIx70O1KtXT2cqG/n1kvy+nS87C9vUgxRJsOYrKuraOwgKynGr7bgGaPQbKG2pYtEi4mYKOm4JnxTSvGHz1zCS1+81NPZMZwsZiNRRsOw0i+P17Vxprkr4CHayTHO9xxfKZjt3TZie82fzEi0jIv0y6PnWum2OyTzRUxIox3U1SilsgBc/9b6uqHW+hGtdbHWujg9fWLk78ZGmXhsfTFRRgN3P7WTQ2dauHfjLnInxbJmSR5vHqrxdJfr6Lbx0q4qrpk7mTQ/mhgISImLItpkGBDU/b8dp+mxOwOAzeUXdtEa27vZX+3cRfNVC9ff5vI6lhWmej0zbDYauHR6KpvL6vuk2JZWNVOUES8fVINk3bICqps6ee+oz7eXUeEZNp4xeB2QyWigMC1uTHbArGjo8HQhHMqFtuih+2B4rqVLdumEX4oyE8ZU87AEi2lY6Zeby5zHpFVF/nW9dHM33fJ3py4zIXpc7NTt88wklJOkYuIZ7aDuT8B619frgddG+fHHvNxJsTy8bjGnz3fwqd98QLfdwaN3FXPvZdOwOzTP76gE4E97z9Bqtfk1s0Y4KaXIcXXAdLM7NM9tq2R5YSq5k2LY3Cv19R+u4vSUuChPJ7XBVJ3v4ERd+6BnVFcVpVPd1MmJeueOoNbOVveSKhI8V87OJCMhmo1hbpjiDtJmZA4+LBmcH0Ddw8nHkspG/4M6d4fMUKZf1rRYx8TOixCBcgZ1ge/UbSmvozAtLuAA1Z1+2eSjxrWj+0JNHThPytS1WsOe4eDLsdo2r93BzzR1su1E75mETUyKNZM7afAB7UKMRyEL6pRSm4CtwEylVJVS6nPAT4CrlFLlwFWu70U/l0xJ4b9unEuM2civ1yxkekY8BalxXDYjnU3bK+mxO9hYUsGsyQksLpCpEIGYlZXAu0drOXKuBYD3jtZS3dTJuuUFrJ6RztbjDfTYnfVzm8vqSIoxc9fyAg6caaZhkPbQgCfwWz3IGVV3WuYWV/D4+AcnOd/Rw9LClBE/N+FkNhpYsySf98vqPN0bw6G8ppVok8GverSijHhON3Z6ZkeNBa1dPTS2d5Of4l8HuSiTgdS4qJCe7a9t6ZLZSCIiJVjMAe/UWW12Sk40BrxLBxd26pp87dRZ++7UZSfH4NBQfX50ZrkG4lR9Ozf97kNueXirZ1wBQEOblVse3srtj5bw9uEawN1JOjmsjXGECJeQBXVa6zVa6yyttVlrnau1flxr3aC1/qjWusj1b//umMLl1uI89v7wKi6feaGXzJ3LCqhpsfKzN49y8EwLdy4rkDeuAP3w2ouJjzbx+ad20tjezcaSCjISorlqdiari9Jos9rYU9nkKk6vZ+X0NK6YmYHWzrbSg9lSXkdWksXTzdCbvJRYpqTGsrm8nneP1vLAXw5zzZzJ3LAgJ9hPdUJbsyQfg1I8uz18u3XltW1MS4/3q4mRezfveG37ELccPe7Ol1P83KkDZwOTUNXlOBya2lYrmTKjTkSgBIsp4EYpu06dp7PHPqwW8kmunboWH0Gds6buQlA3x1XTXVod2KzcUGvp6uHzT+9EKedooLuf3klNSxfdNgdffGY39W1WpqXH87Xn91Ja1URZTWtAXUKFGE/GbKMUwYBOVx+ZlUFOcgy/f/8E8dEmblgogUCgJidZeOSuYmpbrax/Yjvvl9Vx+5J8zEYDy6elYTQotpTXcay2jXMtXawqSmNOThLJseZBUzBtdgcflNezqihtyEB7VVE6Hx6v56vP7WHm5ER+fut8ab0cZJOTLFx1USYv7qzy2dI71Mpr2ijyI/USYLqr7m4spWC6g7r8gIK60NXlNLR3Y3NoSb8UESlxGDt175fXYTYqlhWmBvx4CdEmjAbls1GKs6buQvrlzMkJRBkNlAY4KzeU7A7N15/fy6n6dn57xyIeW19Mm9XGPU/v5Puv7Gf7qUb+5+Z5PP1PS7CYjdzx6DYcGuZKOYOYoCSoiyBGg2LNEmenyxsX5hDfKx9e+G9BXjI/vXke+6ubMagLv9OkGDML8pLZXFbH+670yFUz0jEaFCump7GlvM7rDEGA0upmWrpsfnUoWz0jna4eB1EmA4/etbjPgVUEz53LCmhs7+aNA2dH/bHbrTaqmzopGmTXtreC1FjMRjWmmqVUNDp3DQtS/R/gOznJErKaOncDFkm/FJFoODV1W8rqWVwwqU/tm7+UUiRaTF4bpWit6ei2Exd9YacuymTgouzEPumN4fbg2+W8c6SWDdddzKXT0rgoK5Ff3LaAfVXNvLirii9dPo3rF+SQnRzDI3ctxuoaPSQ7dWKikk+TEWbt0gIOVLdwz+rCcC8lol2/IIemjh46uu1kJV0oqF5VlMav3i4HYFp6nGcMwmVF6fy59CxlNW3MnDywm+HmsjqUgpXTh659WDk9jU/MncznVxX6VW8lhufSaakUpsXxTEklNy7MHdXH9nS+zBy886Wb2Whgaloc5TVjJ6irbOggNS4qoJNHGQkWGtqt9NgdmAOYqeWPC4PHJf1SRB5nTZ3/QV1dq5VDZ1v456tnDvsxk2OjvNbUWW0O7A494ITi/NwkXtpVhcOhx0T2yMt7qrhiZnqfua9XXzyZB26cy7HaNr79sQu/m0X5k/jV7QvYXF7vGa8ixEQjO3URJiUuiofXLR5TrZoj1fpLp/DFy6f1uWxVUTpaOweC9951W+kqVPc12mBLeT3zcpKYFOd7HplbTJSR396xmEX50uQmlAwGxdql+eyqOM+hMy2j+tjlnnEG/u3UOW87tjpgVjR0BJR6Cc6aOq2dH0iDzb0DKOmXIhIlWEy0WW1+d5f84JjzWLM6wPl0vSXGmL12v+zsdqakx/YbozM3J4n2bjsn6sN/cul8ezenGztZ6iX1dO3SfH74qdkDAs9r5mbx40/PHa0lCjHmSFAnRC/zc5NItDjPXq6ecWHXLTs5hukZ8bxfVofdofv819TRzd7TTQEPhxWhd8viPCxmAxtLKjz/vxwjaNntK/22v/LaVqKMBvIDOPlSlBlPZWNHnxpAfx8vFCobOygI8OTR5CT3WIPg19XVtHR5miUIEWkSXMeVNqtzt26ov+0tZfWkxEVxcXbisB8zOcbstVFKu6vLblz/nbo8Zy3avtPhr6srdc2IlXlzQvhP0i+F6MVkNLBiehpvH64dUJy+uiidJ/5xkmnf+4vXnx1OhzIRWkmxZj41L5tN2yvZtN054zElLop3vnUZybFD76r21tjezdW/3Mw3rpzB2qX5g972WE0bhelxA5odDaYoIwGt4XhdG7OzEvneKwc4UdfG8/csG/Uut1abnTPNnRSkBpa26q53C1VQlxYfHfS0TiFGQ6LF2Y2ytauHc81dfPYP27n38mnctXzKgNtqrdns6r48kjTI5FgzpxoGdtTtcO/URffdqZuWHk9slJHSqiZuWjy6Kev9lbpq++bkSFAnhL8kqBOin+9cM4vbLskbUG/whcsKSY2P8po+MynWTLHMDByTvn31TKakxWF3aBrbu3nyw1N8cKyea+dlB3Q/m8vqqGu18oPXDjAlLZZLp/munyyrbWVBXmCvB3enzGO1bbx3tM4ThH54vIEVftRqBtPpxk60xu/B427u1MhQNEupaemSejoRsdw7dZUNHfzLS6Wcae5iw58OembQ9nb4bCv1bdZhzafrLSnG7LVRSrtrt7B/+qXRoJiTk+TZJQunfVXNFKbHeYJhIcTQJKgTop+C1DivHf8yEy3cd8X0MKxIjETv/282u4OXd1exuaxuWEHdpFgzqfHRfOnZ3bx23wqvr5OObhtV5zu5ZXFeQPc/JTUOo0Hx1Ien2HO6iU/OzeLD4/U8U1Ix6kFdpafzZWBBXUpsFCaDCslO3bkWK9lSTyciVIIrOPnGC3s539HDc3cv5Uf/d4gvP7ebV+9bwbT0C/W37trtkWZ/JLuCuv6NTy7U1A38CDg/N4mnt1aEpNlRIEqrmkb9fU+ISCd5LEKICcNkNLCyKI0t5fUB1au506FWFaXz2F3FaA2ff2qn17lTJ+ra0TqwJingbCk+JTWW3ZVNzM1J4ue3zufW4jzePFTDuebQzH7zxTOjLsX/cQbgbE6TkRCaWXW1LV1kSlAnIpR7p66mxcpPb57HpdPSeGx9MVFGA3c/tZPmXvPktpTXMzMzgcwRdnFMjDGjNQO6bra7grr+NXUA83KTsdocHD0XvqZNNS1d1LZapZ5OiABJUCeEmFBWFaVztrnLM3bAH73ToaakxfG7OxZxor6drz+/d0A6rruDpb+Dx3ubn5tMZmI0j6wrxmI2snZpPnaH5vkdlQHf10hUNHQQG2UkLT6wukOA3JRYDp1pCWqTl26bg4b2bjJlRp2IUJOTLBgUntlqALmTYnl43WJOn+/gy5t2Y7M76Oy2s/1UY59GXcPlrhvun4LZ4WqUEtMv/RIuNCYJ5xBy96w8CeqECIwEdUKICcVdp7K5vN7vn3GnQ7k7nF46PY1//9Rs3j5Sy0//drTPbctr2jAbVUBDu90e+PRc3vzGZZ7aNHe9zabtlfTYHQHf33BVNnaQnxI7rAYt183P5si5VvYEcYhxbavMqBORLTPRwrbvXTlg7twlU1L4zxvmsKW8ngf+coRtJxvotjmC0k05OcaZ8tnU2XesQbvVtVMXPTCoy0+JJTnWTGlV+IaQl1Y1YzQoZmdJUCdEICSoE0JMKLmTYilMj/M5c9CbzeV1zMxM6DMjbd2yAu5Yms/D7x/nlT1VnsvLatqYmhY3rHoUi9lIUkzfxgB3LiugpsXK24drAr6/4apoaGfKMIJSgBsW5hAXZeSZkoqgrcfdeEXSL0UkS0+I9nqi5LZL8vnsiik88Y+T/NefDxNtMrBkasqIHy8p1vle4munzltNnVKKuTlJ4d2pq2piRmaC151EIYRvEtQJISac1UXplJxo6DMTzq25o4dfv11OQ5szkOjstrPj5PkBneiUUmy47mKWFabwry/t51//WMp3Xy5lZ0UjRRkJQVvrR2ZlkJMcw5MfnhqVuXV2h+Z0Y2fATVLc4qNNfHpRLq+XnqWxfeDg4+FwN16R9EsxXn3/ExexqiiN8to2lkxNwWIeeUDj2anr6B/UeR8+7jY/N5mjNa1e3x9DTWvN/upm5kvqpRABk6BOCDHhrCpKo6vHwa6K830ut9kd3Pfcbn7+VhlffGY33TaHMx3K7mCVl050ZqOB396xmPm5Sbx7tJa3D9cSbTJw1ezMoK3VaFB8buVUSk408uiWE0G7X18qGtrptjsoTB/eTh04dxe7bQ5e3Hk6KGv64Fg90SYD+cMMNIUY60xGA79Zs4jlhancdklgnXN9SfKkX/YN6tq7bUQZDT6zCZYWpmB3aN46NHrZAW6nGztp6uhhrgR1QgRMRhoIISacZYWpmI2KzWV1fdpm/+efD/PBsXo+vTCHl/dU88PXDhAbZSLKZGCpj3SolLgoXrz30pCu97MrprCr4jw/fuMIRRkJXDErI2SP5U67mpebPOz7mDk5gSVTUnh2WyV3ryoc0QDl1q4eXt1TzXXzs4mPlkOWGL+SYs1sumdZ0O4v0RXUtfQL6jq77QMGj/e2YloaeSkxPFNSwafmBzb6ZaT2uWr55o/g/UeIiUp26oQQE05ctInFBZP6NEvZtL2SJz88xedWTuV/b1vAfVdM4/kdp3l2WwVLg5QONVxKKX52y3xmZyXy1U17OFYbunbjpVXNWMyGgEcy9Hfn8gIqGzvYHEDtojev7Kmmo9vOncsKRnQ/Qkw0FrORGLORpo6BjVK8jTNwMxgUdywtYNvJRspqRne0QWlVE1EmAzMyg5fCLsREIac9hRAT0qqidH76t6MsfeDvANS3dbN6RjrfvWYWAN+6aiZlNW28dahmQD1dOMREGXn0rmKu+80/+NxTO3ntvhWeluXBVFrVxMXZSZhGOHj44xdPJi0+ioffP86qonSMPnbr/ririjcPnuP36xYPaCKhtWbj1grm5SYxP0/O3AsRqCTXAPLeOrptPuvp3G4tzuN/3yrj2ZIK7r9+jt+P19rVw5pHS6hrddYkmwwGvnHVDG5enOvzZ76yaQ/bTzYAzvq/i7ISiTLJnoMQgZK/GiHEhHRLcS53LS/gipkZXDEzg8+vnMqv1yz0BDMGg+IXty3gqx+Zzs2Lg1PjMlLZyTH8ft1izjZ1cd9zu4M+5sBmd3DgTHNQ5kNFmZwf5kpONA4Y++Bmd2h+8VYZbx6q4WR9+4Drt59spLy2jTuXyi6dEMORHGse0Cilvds+ZFCXEhfFJ+dm8dLuatqttkFv29ure6o5UN3C8sJUrpiZQWp8FN95qZStxxu83v5AdTP/t++MM618ZgY3LszhW1fN8PvxhBAXyE6dEGJCykiw8KMhzkDHR5v45sdmDnqb0ba4YBIPfHou335xH//5+qGAzqIPpby2ja4eR9DqWdYuyefQmRYefv84MyfHc+PCvmfr3ztaS3VTJwCby+ooTO+b8rmxpIJEi2nU63qEGC8SY8wDGqV0dtu8jjPo785lBbyyp5rX9p5h7dL8IW+vtWZjSQVzc5L4xW0LUErR0tXDjQ/9gy89u4vX7ls5oNnRs9sqsJgNPHTHogHjXIQQgZGdOiGEiDA3L87l7lVTeWprBc9tqxzy9r5GC/RPy9rvaZISnM5z/cc+7Kns2230mZIKMhKiyUuJYUu/YfC1rV389cA5binOk3lVQgxTcox5QKOUdqvd6+Dx/hblJ3NRViIbSyq8jlPp/76y49R5ymrauHNZvieVOtFi5rH1l+DQcPfTO2nrtevX0tXDq3vOcP38HAnohAgCCeqEECICfeeai7h8Zjo/fO0AJSe8pzYBPLL5OIv+4y2e3943+PvTvjMs/NGbvF56xnPZvqomEqJNwx487o177ENmYjRf2LiLc83OmXOnGzt4r6yO25fkc/mMDLaeaKDbdiGd9IUdp7E5NHf4sUMghPAuKWZg+mVHt40YP3bqlFKsW1bA4bMtvLy7us91f9xVxaL/eIvfvnfMc9nGkgoSLCaum5/T57ZT0+J4aO0ijtW18fXn9+JwOAPEl3dV0dljZ91ySa8WIhgkqBNCiAhkNCgeXLOQgtRYvvjMLk43dgy4zTtHavjxG0eIjzbxg9cOsP1kI+BshvLPL+7DoeHxD056bl9a1czc3KQRjSDwJiUuisfXX0K71cY9G3fS1WPn2W2VGJRizZI8VhWl0dFt98wNtDs0z22rZOX0tAEpmUII/yXHDmyU0t5tJ87P3e9binNZXpjKd1+5sNO+q6KR7728n/hoEz/921HeOlRDXauVvx44yy2Lve+sryxK4wefvIi/H67hZ28eRWvNM9sqmZ+XzJwcmUknRDBIUCeEEBGqd2rT55/qm9pUXtPKVzftZXZWIm9/6zLyJsVy7zO72FXRyN1P7yQtPpovXzGdPZVNHKhuxmqzc+Rcy4jm0w1mRmYCv7p9Ifurm/nWi/t4Yedprrwog6ykGJZPS8VkUGxxjT9450gtZ5q7uHOZ7NIJMRLJsVF09tjp6rF7LuvstvtVUwfunfZFTE60cM/GXeyqOM8XNu4iK9nC3795GfNykvj683v48RuH6bFr7hjkb3b9pVNYsySP3753nB++dpBjtW2sk1ElQgSNBHVCCBHBeqc23fP0Tn7zTjm/eaeczz+9E4vZOQYhM9HCo+uL6bE7uOl3W2ntsvHY+mLuXl2IxWzg2W0VHDnbSo9dMz9I9XTeXDk7k3++eiZ/Lj1LY3s365ZNASDBYmZR/iTPTLuNJRVkJkZz5UWZIVuLEBNB/wHkWmvau21+1dS5TYqL4rH1xXR227npdx9i7XHw+PpiJidZ+P26YuKiTby8u5oV01OZNsjOulKK+6+bw5KpKWwsqSApxsy187JG9gSFEB4S1AkhRIRbWZTG/dddzI5TjfzszTJ+9mYZTR09/H7dYrKTYwCYlh7PQ2sXkZ4QzS9uW8BFWYkkxZi5fn4Or+45wwfHnI1K5oV4HtwXL5vG2qX5LJmSwqXTUj2XrypK40B1C7srz7O5rI41S/JHPCtPiIku2RXUuVMwu3ocaE3AzYdmZCbw4JoFpCdE8+CahUzPcA4Hn5xk4ZG7islOsvDFy6YPeT9RJgO/u2MRF2Ulcu9l07CYpQmSEMGivHU0GmuKi4v1zp07w70MIYQY02x2B+53dINSXgd+a637DPk+UN3Mtb/+gPhoE9EmAzv/7coBQ8BDof869p1u4vqH/sHUtDgqGzv48DsfITPREvJ1CDGebS6r464ntvPivcu5ZEoKDW1WFv/n37n/uotZf+mUgO+v/9/tUJcHej9CiMEppXZprYu9XSenQYUQYpwwGQ2YXf95C+iAAR+k5uQkMT8vmTarjXm5SaP2QcvbOpJjzZysb+djszMloBMiCJJjXTt1rg6YHd3O2rqhho/74uv9IdD3DQnohAg+CeqEEGKCczcrCFWTFH8YDYoV09P6rEcIMTLJMVEAngHk7d3OZkr+NkoRQkSOsPxVK6U+DvwKMAKPaa1/Eo51CCGEgGvnZbG78jw3LMwZ+sYh9LmVU5mcaGF5r1o7IcTwuYd6N3U4B4W3W107dQE0ShFCRIZRD+qUUkbgIeAqoArYoZT6k9b60GivRQghBFjMRh64cW64l8Gi/Eksyp8U7mUIMW4kWEwodaH7Zacr/TJOduqEGHfCkX65BDimtT6hte4GngeuD8M6hBBCCCHGLYNBkWgxe0m/lJ06IcabcAR1OcDpXt9XuS4TQgghhBBBlBxr9ow06JCgTohxy7hhw4ZRfcD777//YmD6hg0b/uT6fj6Qs2HDhjd6304pdc/999//+/vvv/8ek8mUvXDhQmJiYjh48CDHjh0jNTWVrVu30tHRQXNzM6WlpcTHx7N3715OnTpFUlISJSUlWK1WGhoaOHDgAMnJyezYsYPq6mpiYmLYvn07NpuNc+fOcfDgQc991tTUYDabcY9ROCfQVM0AAAmrSURBVH36NIcPH/ZcX1/vnOe0a9cujEYjJ0+e5MiRI57rm5qa6OnpYc+ePURFRXH06FHKyso817e2ttLZ2cnevXvlOclzkuckz0mekzwneU7ynEL2nA7XdxPfVk18z3nOdCiaKw6xKDeB9pbzEfucxuP/J3lO8pz8eU5PPfXU2Q0bNjziLcYa9Tl1SqnlwAat9dWu778LoLX+sa+fkTl1QgghhBCBW/f4Nlq6bLx23woe3XyC//rLYfZv+BgJFnO4lyaECNBYm1O3AyhSSk1VSkUBtwN/CsM6hBBCCCHGteTYKM42dfLy7ip2VZwHZKSBEOPRqP9Va61tSqkvA3/DOdLgCa31wdFehxBCCCHEeFeQEsv/7TvDN1/YB8DkRAtGgwz/FmK8GfX0y+GQ9EshhBBCiMA5HJqq851onJ/3UuKiJPVSiAg1WPql7L8LIYQQQoxTBoMiPzU23MsQQoRYOGrqhBBCCCGEEEIEiQR1QgghhBBCCBHBJKgTQgghhBBCiAgmQZ0QQgghhBBCRDAJ6oQQQgghhBAigkXESAOlVB1QEe51eJEG1Id7EWJck9eYCDV5jYlQkteXCDV5jYlQG0uvsQKtdbq3KyIiqBurlFI7fc2KECIY5DUmQk1eYyKU5PUlQk1eYyLUIuU1JumXQgghhBBCCBHBJKgTQgghhBBCiAgmQd3IPBLuBYhxT15jItTkNSZCSV5fItTkNSZCLSJeY1JTJ4QQQgghhBARTHbqhBBCCCGEECKCSVAnhBBCCCGEEBFMgrphUEo9oZSqVUodCPdaxPijlMpTSr2rlDqslDqolPpauNckxhellEUptV0ptc/1Grs/3GsS45NSyqiU2qOUej3caxHjj1LqlFJqv1Jqr1JqZ7jXI8YfpVSyUuqPSqkjrs9ly8O9Jl+kpm4YlFKrgTbgaa31nHCvR4wvSqksIEtrvVsplQDsAm7QWh8K89LEOKGUUkCc1rpNKWUGPgC+prUuCfPSxDijlPomUAwkaq2vDfd6xPiilDoFFGutx8pgaDHOKKWeArZorR9TSkUBsVrrpnCvyxvZqRsGrfVmoDHc6xDjk9b6rNZ6t+vrVuAwkBPeVYnxRDu1ub41u/6TM3wiqJRSucAngcfCvRYhhAiUUioRWA08DqC17h6rAR1IUCfEmKaUmgIsBLaFdyVivHGlxe0FaoG3tNbyGhPB9kvgXwBHuBcixi0NvKmU2qWUuifcixHjTiFQB/zBlUb+mFIqLtyL8kWCOiHGKKVUPPAS8HWtdUu41yPGF621XWu9AMgFliilJJVcBI1S6lqgVmu9K9xrEePaCq31IuAa4D5XeYwQwWICFgG/01ovBNqB74R3Sb5JUCfEGOSqc3oJeFZr/XK41yPGL1cqyXvAx8O8FDG+rACuc9U8PQ98RCn1THiXJMYbrfUZ17+1wCvAkvCuSIwzVUBVr0yWP+IM8sYkCeqEGGNcTSweBw5rrf833OsR449SKl0plez6Oga4EjgS3lWJ8URr/V2tda7WegpwO/CO1vrOMC9LjCNKqThXMzFcKXEfA6QruQgarfU54LRSaqbroo8CY7ZpnSncC4hESqlNwOVAmlKqCvh3rfXj4V2VGEdWAOuA/a6aJ4Dvaa3/EsY1ifElC3hKKWXEeXLvBa21tJwXQkSSTOAV53lQTMBzWuu/hndJYhz6CvCsq/PlCeCzYV6PTzLSQAghhBBCCCEimKRfCiGEEEIIIUQEk6BOCCGEEEIIISKYBHVCCCGEEEIIEcEkqBNCCCGEEEKICCZBnRBCCCGEEEKEiFLqCaVUrVJqyLEbSqlfKKX2uv4rU0o1+fMYEtQJIYQYk5RSTyqlbvbztn4fMIOwrg9d/05RSq3tdXmxUurBEDzeKaVUWrDvVwghxKh5Evi4PzfUWn9Da71Aa70A+DXwsj8/J0GdEEKIiOWatQcBHDBHSmt9qevLKcDaXpfv1Fp/dTTWIIQQInJorTcDjb0vU0pNU0r9VSm1Sym1RSk1y8uPrgE2+fMYEtQJIYQIO6XUXUqpUqXUPqXUxl5XrVZKfaiUOuHetVNKXa6Uelcp9RywH7wfML08xpNKqYddB88ypdS1rsstSqk/KKX2K6X2KKWucF1+sVJquysFplQpVeS6vM11lz8BVrmu/4ZrXa+7bpOilHrV9XMlSql5rss3uHYV33M9J3+DwH92rWW7Umr6YM9HCCFERHgE+IrWejHwbeC3va9UShUAU4F3/LkzU9CXJ4QQQgRAKXUx8H1ghda6XimV0uvqLGAlMAv4E/BH1+VLgDla65MBPtwU4DJgGvCuK0C6D0BrPdd1pvRNpdQM4F7gV1rrZ5VSUYCx3319B/i21todHF7e67r7gT1a6xuUUh8BngYWuK6bBVwBJABHlVK/01r3DLHuFq31EqXUXcAvAXcAN+D5aK27/P91CCGEGG1KqXjgUuBFpZT74uh+N7sd+KPW2u7PfUpQJ4QQItw+gvPAVQ+gte694/aq1toBHFJKZfa6fPswAjqAF1z3V66UOoEzwFqJs24BrfURpVQFMAPYCnxfKZULvKy1Lg/gcVYCN7nu8x2lVKpSKsl13Z+11lbAqpSqBTKBqiHub1Ovf38xxPPZG8A6hRBCjD4D0OSqm/PldlwnHf29QyGEECKcFKB9XGftdzu39mE+Vv/H0f3u98IVWj8HXAd0An9z7bj5y9t9uh+793Oy498JVu3H196+F0IIMcZorVuAk0qpWwCU03z39UqpmcAknCcX/SJBnRBCiHB7G7hVKZUKznq0ED7WLUopg1JqGlAIHAU2A3e4HnsGkI8zLbIQOKG1fhBn6ue8fvfVijOF0pve93k5UO86iPuklHpbKZXj4+rbev3b+yDv7fkIIYQYQ5RSm3C+d89USlUppT6H8xjxOaXUPuAgcH2vH1kDPK+19vtEnaRfCiGECCut9UGl1H8B7yul7MAe4DOB3IfrgHk5kKaUqgL+XWv9uJebHgXex5nyeK/Wuksp9VvgYaXUfsAGfEZrbVVK3QbcqZTqAc4BP+p3X6WAzXVAftK1brcNwB+UUqVAB7B+iPUbgOn4bvYSrZTahvNk7JrBns9gjyOEEGL0aa3X+LjKa9dmrfWGQB9DBRAACiGEEBFLKfUk8LrW+o9D3Xa0KaXmAP+ktf5mAD/zJGP0+QghhBhdslMnhBBChJnW+gDgd0AnhBBC9CY7dUIIIYQQQggRwaRRihBCCCGEEEJEMAnqhBBCCCGEECKCSVAnhBBCCCGEEBFMgjohhBBCCCGEiGAS1AkhhBBCCCFEBPv/94tvhCDZ2wsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x216 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the digitized E1 to confirm that digitization was successful.\n",
    "\n",
    "plt.figure(\n",
    "    figsize=(15,3)\n",
    ")\n",
    "\n",
    "loc_eig = bioframe.slice_bedframe(digitized, 'chr1:10M-60M')\n",
    "plt.plot(\n",
    "    loc_eig['start'],\n",
    "    loc_eig['E1.d']\n",
    ")\n",
    "plt.axhline(0,ls='--',lw=0.5,color='gray')\n",
    "plt.ylabel('E1, digitized')\n",
    "plt.xlabel('chr1 position, bp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:57:24.250670Z",
     "start_time": "2018-09-25T16:57:03.985462Z"
    }
   },
   "outputs": [],
   "source": [
    "# Calculate the decay of contact frequency with distance (i.e. \"expected\")\n",
    "# for each chromosome.\n",
    "\n",
    "expected = cooltools.expected.cis_expected(c, regions, use_dask=True)\n",
    "\n",
    "# Make a function that returns observed/expected dense matrix of an arbitrary\n",
    "# region of the Hi-C map.\n",
    "getmatrix = cooltools.saddle.make_cis_obsexp_fetcher(c, (expected, 'balanced.avg'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:58:45.310373Z",
     "start_time": "2018-09-25T16:57:24.253749Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/s1529682/Projects/cooltools/cooltools/saddle.py:144: RuntimeWarning: invalid value encountered in true_divide\n",
      "  return obs_mat/exp_mat\n",
      "/home/s1529682/Projects/cooltools/cooltools/saddle.py:144: RuntimeWarning: divide by zero encountered in true_divide\n",
      "  return obs_mat/exp_mat\n"
     ]
    }
   ],
   "source": [
    "# Compute the saddle plot, i.e. the average observed/expected between genomic \n",
    "# ins as a function of their digitized E1.\n",
    "\n",
    "S, C = cooltools.saddle.make_saddle(\n",
    "    getmatrix,\n",
    "    group_E1_bounds,\n",
    "    (digitized, 'E1' + '.d'),\n",
    "    contact_type='cis')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-09-25T16:58:45.511782Z",
     "start_time": "2018-09-25T16:58:45.314203Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.colorbar.Colorbar at 0x7fdb7cdb6d68>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAD8CAYAAAARze3ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2de7RkV13nP9963Uc/0ul0EjoPCcsVXCILgvYEXLjGgIRpMmDEASbBYVgO2OAiDqjLRUZmBOO4jIoiKhIbiASXElFAMhgJIcoEZIAEDHkAISFE6KTzTvp969bjN3+cc0NVnd++dapv3eq+dX+ftc66Vfvss88+p+r+ap/93b/fT2ZGEARBkKZyrDsQBEFwvBOGMgiCYAhhKIMgCIYQhjIIgmAIYSiDIAiGEIYyCIJgCGEogyBYNSRdKekhSbcn9kvSH0u6W9Ktkn60Z99OSXfm+y6dXK+LhKEMgmA1+SCwc5n9LwHOzrddwHsBJFWB9+T7nwFcLOkZq9rTZQhDGQTBqmFmNwKPLVPlQuBDlvFFYIuk7cC5wN1mdo+ZLQJX53WPCbWVHCxpJ/BuoAq838wuX67+CaraKdRXcspVpdKQX16v+uXVYv1k3dpo5ao4v2E1/+NS1W/Dqn59qxTrdyt+3a788rYV2+iY/7vb6fj3tdl2i2m3uoWyVqvj1m01W255t+PXn3aahx+gtbjPv+El+bHKBttv5e7f3TTvABZ6inab2e4RTnc68L2e93vyMq/8uSO0O1aO2lD2DI3PJ7uImyRdY2ZfTx1zCnXeVX3q0Z5y1Znb3nDLN52xwS2f3TxTKNtw8ia37vy2zX4bJ5/oltc2Fs9ZO9Gvq80nuOWdE09xy1vzWwplC3N+2wcaJ7nlj7WKbexrzrl1Hz/k/zh+9wG3mIcfbhbKHrjvCbfuA/fsdcsP7TvgNz7lfO1zu1bcxn7r8Ee1cv+nL21/a8HMdqzgdJ5Rt2XKjwkrGVE+OTQGkLQ0NE4ayiAI1gAC1UsOShNPBSOwBziz5/0ZwP1AI1F+TFiJoSw1NJa0i2ySlpNX9qQfBMEEUEVU5/zpnAJHVny6a4BL8oHWc4F9ZrZX0sPA2ZKeBtwHXAS8esVnO0pWYrlKDY3z+YrdAGdrNkIVBcHxjqBSW9E05/ebkj4MnAdsk7QHeDtkQoWZXQFcC1wA3A0cBn4+39eWdAlwHZkGcqWZ3TGWTh0FKzGUqSFzEARrmVEevYdgZhcP2W/AmxL7riUzpMeclRjKmziOhsajUN/sP1bMnzLrltdm/PqecFOb9YWL2ga/7eqcX16ZKQpLmvPFEpv3BaTOjC9CNWeL4s+Rui82tcwXuA61ikLWYwf9az+04P/TpZTshYWikn3w8UN+2+tUtFlNJI1tRDktHLWhPN6GxkEQjIkxjiinhRWpK8fT0DgIgjExxjnKaSFk6CAI+pCg2ginvV7CUAZBMIBQJUaUvaxLQ5maf/FcEmE5N8Ni/WrDv6W1lGhTT7glzhTFEuq+sIISv/6JcpPjwuiUAbS6if6puNKrmuhGqjyFFT0Yacz5196Yde4TsLhQ9O4JSiLQqB/alLMuDWUQBGlEetCwXglDGQRBP/KfltYzYSiDIOhDUog5A4ShDIKggBvmbx0ThjIIgn7i0bvA1BvKxknFS5zd5iuo9XnfBa8+Vz7YcEoh77Z9dz3rJuKEmFPeTQRT7fjBa6sLvntfo150hVys+e6R84mRxeFqUcXfPOf3o9P1718j8Xh30rZiXxaOLLp122f4MTcf+PYet9y8+xoMoBBzBph6QxkEwWgoRpQFwlAGQVAg5ij7CUMZBEE/EtV6GMpewlAGQdBHPHoXCUMZBEGBePTuZ2oMpaduA9Q3FVXomc2jqdu1mYS/s/NlSqrYXceBGbAR0qpay1eVGTE1a6VTVJBrXV9VblZ8NbxRKWaVqlf9+5dyG9680R+1tNvFA2YTvt6pdLWzm/ygxUf2H/Q7E3yfGFEWiJ+NIAgGyKIHldlKtSbtlHSnpLslXers/zVJt+Tb7ZI6krbm++6VdFu+7+YxX2hppmZEGQTBeJDS64FHb0tV4D3A+WR5tm6SdI2ZPZnW2sx+H/j9vP7LgF82s8d6mnmBmT0ylg4dJWEogyAoMMYF5+cCd5vZPQB5WtoLga8n6l8MfHhcJx8X8egdBEE/Guuj9+nA93re78nLnNNqHtgJfLSn2IBPS/qKpF1HeUUrJkaUQRAUGEH13jYwd7jbzHb3NuUck/IjfRnwLwOP3c83s/slnQJcL+mbZnZj2c6NixUZSkn3AgeADtA2sx3j6NRyjKJuAzQ2lffTTpHy067Uil+m1K9sd0Rl2lPDK9XEx9VOqOGpyOcOVUcJB6hVE37kKvZvY33Brbs45382+w8nfMCdCPRbtvoR4ptnbnPLH3RLQ/Uuw4jrKB8Z8n+/Bziz5/0ZwP2Juhcx8NhtZvfnfx+S9HGyR/mJG8pxPHq/wMzOmYSRDIJgMozx0fsm4GxJT5PUIDOG1xTOJ50A/CTwiZ6yDZI2Lb0GXgzcPobLG5l49A6CoB9pbKq3mbUlXQJcB1SBK83sDklvzPdfkVd9OfBpMzvUc/ipwMclQWar/trMPjWWjo3ISg3l0kSrAX8+MDcBQD4Buwvg5LDLQbAG0Fg9c8zsWuDagbIrBt5/EPjgQNk9wLPH1pEVsFLLNXSiNTeeuwHO1mwEAwyCtYDCM6eXFRnKcU201jcXh/mplLKpoLspt0S/ri8MpOZc6vN+SlTv8aQ259etzvj9Trkwqlb8aKyVSME657vrqe0LNHLywVpC+KmY379GpSjyLCZS21ac1LYAicy+zM46aYAP+J9NNeEfObfRd72sOemB262iO+Z6JoJiFDnq8fXxNNEaBMF4UaVSalsvrGREedxMtAZBMEZU3o97vXDUhvJ4mmgNgmC8jEv1nhZChg6CoI+YoywShjIIggEE62j+sQwTNZSVhpjbXlR/508pqtCp6CXjSCmb+rWcP2mjWz6Ox5DkxHcqoG+7qMQqpZA3D/ttO8o5QG2h6MY3kwi6a/KvveOox6lvU0oN37opobRXivWbi34/mk1f3W4n3E7bT91eKNt79/ecmusbxfKgPmJEGQRBP4pUEIOEoQyCoB8JhZjTRxjKIAgKhJjTTxjKIAj6EEIjhOhbD4ShDIKgHwExouxjsqp3vcqmM4p+ybWZ4nxISmkeNaWsF3R3FN/t5cqrjrNyMl1tglR9LzWtJZRIzfu+3rR8X29mnKDAHT9Ab6N1yC3vOMr0YsW/rzX5yv7Ghn/ORSdd7Qkb/RHOwoL/fWgt+mp4q1lcTTC32V/tsJ6D/IaY00+MKIMgKBBzlP2EoQyCoB8JVUP17iUMZRAEReLRu48wlEEQ9CEpPHMGCEMZBEGRGFH2MVnVuypmNxeV0Q0nbyqUjTqZnFLpvJSyo/pue+o2+NHMk5HMEz7d1XlfncXz607NG6WU9lTUckcNry4ecetaxT/njOMD3pr1Ve/Zqh9RvtX1256tF6+94XyOABs3+OXttt+Xw4eLSvv8Jn/VwPpWvcc3opS0E3g3WXKx95vZ5QP7zyPLvvidvOhjZnZZmWMnRYwogyDoJ4uzNqamVAXeA5xPluP7JknXmNnXB6p+zsxeepTHrjoxvg6CoICq1VJbCc4F7jaze8xsEbgauLBkN1Zy7FgJQxkEQT9LnjllNtgm6eaebddAa6cDvXHs9uRlg/y4pK9J+kdJPzLisatOPHoHQTDASHm9HzGzHcs2VmRwUv2rwFPN7KCkC4C/B84ueexEGGooJV0JvBR4yMyemZdtBf4GOAu4F3iVmT0+rK1KveoKN7XZohtaSkBJCTEpV0BvUjqVUjZFqm1PuEmFpxo1cK8n3Nii75JozQX/nInAvd7VqO6LUJVOIpWr4zk40/EDCG+u+cJAx/x71bXivWrP+/evkxCELNF2uz1fKDu4/wS37uKCnx74wKNPuOVTxfiWB+0Bzux5fwZwf28FM9vf8/paSX8maVuZYydFmZ+NDwI7B8ouBW4ws7OBG/L3QRBMAyJbHlRmG85NwNmSniapAVwEXNN3OukpyhduSjqXzC49WubYSTF0RGlmN0o6a6D4QuC8/PVVwGeBt46xX0EQHDPG58JoZm1JlwDXkS3xudLM7pD0xnz/FcArgF+U1AaOABeZmQHusWPp2Igc7RzlqWa2F8DM9ko6JVUxn9zdBXDahsSawSAIjh/E2JYHQfY4DVw7UHZFz+s/Bf607LHHglVXvc1st5ntMLMdW2f9ebAgCI4nSire6yjC0NGOKB+UtD0fTW4HHhpnp4IgOHYIIsL5AEdrKK8BXgtcnv/9RJmDKrUq89s2FzuxoejiVpvz3d68QLzZDl897jqugCk3w1GDlXpuiak2qgn3vhRjmSPqJtLbtotufJVFXzkn4cJYd8qrHV+Vb8/797tR8etvbhTva6vr39fZhl/emvVHO/vrxfobN/tTQodPLK7QgHWgekeE8wJllgd9mEy42SZpD/B2MgP5EUmvA74LvHI1OxkEwSQZnwvjtFBG9b44seunxtyXIAiOFyJwbx/hmRMEQT9jDIoxLYShDIKgSMxR9hGGMgiCIjGi7GOygXtrVWZPPrFQXnUU7ko94aecSu/qBbpNkKybUM6Tvt5e0N2U73aCSirQr9PHpCrvqNgA5sfiRQ2nbb8qlYRyXnP+kbp1X9mfqfuq8ua6f9ZWtdhOq5FIJbzJ/2waNf++er7h7bYfuLeSGFUdfPyAWz5VanikgugjRpRBEPQjRSqIAcJQBkFQJLF+dr0ShjIIgn5iRFkgDGUQBEVijrKPMJRBEBQJ1buPiRpKVSrUNhYVRk/51UzCN9pGiwTvqsepyN9tP5q3tXxVeZSUsinf7ZQCX2kU70lSrU95UaziWrhK24n+nRiF1Bf9yOedqq9MtyvF8o21RDrdhF7fbPv3ZMvGYv3DC37dVsuPN7DttG1u+fSo3ooR5QAxogyCoA8TWLgw9hGGMgiCAcKFcZAwlEEQFAlD2UfcjSAICphUaiuDpJ2S7pR0t6RCIkJJPyfp1nz7gqRn9+y7V9Jtkm6RdPNKrytPZPbTkl4m6Sllj5vsiLJWo3Zi0YVRc46LWyJ9aioYbUpwqVSLl2gtPw2pEmJJ8gsxQkrZFCm3RE+4qcwVU63CMmJTqrxe3oXRS/eb7SiWp351azVfzJnDd/eszBRFNTX8lLIt5/MF2JpIz1SvFgWaI01fOOx0nJy8wMObfJHnxKcURZ7HH3jE78jxzBijB0mqAu8BzidLP3uTpGvM7Os91b4D/KSZPS7pJcBu4Lk9+19gZiu+kZJeD/wG8E9kX/k/kXSZmV057Nh49A6CoMj4VO9zgbvN7J6sWV1NlsX1SUNpZl/oqf9Fsvzdq8GvAc8xs0fzvpwEfAEYaijj0TsIggGEVaulNrLMBzf3bLsGGjsd+F7P+z15WYrXAf/Y896AT0v6itP2qOwBeiOaHBjoW5IYUQZB0M9o6WofMbMdQ1obxF0MLekFZIbyJ3qKn29m9+cpsa+X9E0zu7Fs5wa4D/iSpE/kfbgQ+LKkXwEwsz9MHRiGMgiCAjY+1XsPcGbP+zOA+wcrSXoW8H7gJUuPxgBmdn/+9yFJHyd7lD9aQ/ntfFtiKSmin0WuhzCUQRAMMFbPnJuAsyU9jWxEdxHw6r6zST8AfAx4jZl9q6d8A1AxswP56xcDl62gL79rZn3pRiVtKyMUlcnCeCXwUuAhM3tmXvYO4BeAh/Nqv25m1w5tq1pFm4vqpc07Bj31i9YZwZ0Q/KC2c36gVjV9dVbzfn2cgL7WTKR9TZEIuusq6gkVOxnQt+artngunClPjETQYu/zUcIFVJ2Ea6iTOjZFteu3MSv/fjfxr13OU9+GWf8aD834xmLLVn/1waH9xe/2wmG/f0f2H3TLjxfGNaI0s7akS4DrgCpwpZndIemN+f4ryJTok4A/U2ag2/nj/KnAx/OyGvDXZvapFXTny5J2mdkXAST9J+B3gKcPO7DMiPKDwJ8CHxoof5eZvXPEjgZBcLwjjTUeZT6Iunag7Iqe168HXu8cdw/w7MHyFfBzwJWSPgucRmacX1jmwDLpam+UdNYKOhcEwRrCWGbt8BrGzG6T9NvAX5Ip3v/ezPaUOXYl4+tL8pX0V0oqriLPkbRraenAw/uO78eNIAhyVCm3rSEkfQB4C/As4OeB/yPpTWWOPdorfS/wg8A5wF7gD1IVzWy3me0wsx0nn7DxKE8XBMEkMVRqW2PcTubl8x0zuw54HvCjZQ48KkNpZg+aWcfMusD7yCT7IAimAmGqlNrWEmb2LmBW0g/l7/eZ2evKHHtUy4MkbTezvfnbl5NZ6uEdrdbonHhKobwz46jKiQ+huuCnCk3iqrMJf+xEQF9aifpe26k2Uj7qiZSybtDdlO92St1O4V1Pwmearn/t7irihHJeSfjtNyzh6931Viq4ValU/fuquq9k1yvF+i0nhS1Ap+vfk30H/Pu9f0tRDd//mK+QH++q91p7rC6DpJcB7wQawNMknQNcZmY/PezYMsuDPgycR+aqtAd4O3BefhID7gXecNS9D4LguMIkutOZhfEdZE+/nwUws1vy9Z1DKaN6X+wUf2CEzgVBsNaYQtWbbH3mPvVfW6ncMuGZEwRBgbU2/1iS2yW9GqhKOhv472TRg4YylXcjCIKVUE7xXoOq9y8BPwI0gb8G9pEtFxpKjCiDICgwTSNKSf8D+JSZ/SvwtnwbiYkaSqtUac1vKZQ3Zx3/b/mTyY26L39WOuUjiyuhttYWEkrkTCKlrKMepyY8lPDpViPho+7hRCYHRlfrPT/ylI96wodec8497Pr3tZq43zbjK8JeSuL5xKqBduL7YPP+P3rX8Ys/wUmXDFCRH/n84El+/SeeKPaledpWt27zkL/c4dC+EVd1rAZi2uYovwO8OU8x8TWyeJefNrPHyzYQI8ogCPowRDcxUFmLmNnVwNUAkp4D7AQ+lqep+AzZaPPLy7URhjIIggLT9OjdS/74/a/A70jaTJbL5/XAsoZyOu9GEAQrYhrFHEmvlLQpf/0/ySKjfcfMhqaYCEMZBMEA0+nCCPyvPAjwTwD/AbiKLG7FUCb66N2t1FiYKwYaOlLfXKybmCNZrPmT97WEq13VEXlSH/BM1XdNqySCBVcXixPySrjrVRZ9sWSU3+Rk3VTQ3YRboivcJAQX6ok2POEmIbgkSdzXSruYTthqifuaEPHqHT8lcb1eFGg21Hxhpd31vyeb5v3vyZYtxfJWyw/6vP8xP/3ukYN+8OhuKjD1KjGNYdaApZv4H4H3mtkn8iDkQ4k5yiAI+jBNl5jTw32S/hx4EfC7kmYo+VS95sbOQRCsPtM4Rwm8iiwlxU4zewLYSpbreyhhKIMgKDDOOUpJOyXdKeluSZc6+yXpj/P9t0r60bLHjnRNZofJgvi8RNIvAdvN7NNljg1DGQRBgXGNKPO1iu8BXgI8A7hY0jMGqr0EODvfdpELLCWPLY2k3yATcE4CtgF/kavfQ4k5yiAI+rBc9R4T5wJ354nCkHQ1cCHw9Z46FwIfMjMDvihpi6TtwFkljh2Fi4HnLKWslXQ58FXgfw87cLKqt2ocaJxUKG9ZUdFsJYKmzidSszYrCTW8WlRWK5YIoptym2wd8us7MfsqidSsqax2lRGUYnnBfCGdUjaxEsB1S0yo28k0wN4/Uitx7Qu+qqxq+f5VE/epklipkHRTnSkq/pXZ4ncSoDLjt3Fkg3/Ok070lPlZt25zwT/n3Aa//pFDxX43nVS4tdTnOCIjzD9uk3Rzz/vdZra75/3pwPd63u8BnjvQhlfn9JLHjsK9ZB/I0o2bAb5d5sAYUQZBUKBbflbukTwHdwo3GH7JOmWOHYqkP8mPawJ3SLo+f38+8PkybYShDIJgAGHjky/2AGf2vD8DuL9knUaJY8uwNOL9CvDxnvLPlm0gDGUQBH0YIz16D+Mm4Ow85cJ9wEXAqwfqXEOW/vpqskfrfWa2V9LDJY4dipldtfRaUgN4ev72TjNLJKLqJwxlEAQFxmUozawt6RKy9YtV4Eozu0PSG/P9VwDXAhcAdwOHyXJuJ4892r5IOo9M9b6X7LH+TEmvNbMbhx0bhjIIggLjXExuZteSGcPesit6XhvwprLHroA/AF5sZncCSHo68GHgx4YdWCYL45nAh4CnAF0yVevdkrYCf0Mm4d8LvGpYIMy2VXmsVQzce6hV9L+V/Dnbw1VfFWxUfMW1qqJa2qj4o+1OQjHsVPzyGU8lT2SOrSdU71piGYbn75wMpppoI/VV94Luur7by7SdVMNHIaHWu0GOE4GP0z7g/veh0dxfKJur+m10Gv7nvrHuq/WbHR/whab/uZ+wxf8eW+qeOCseGrPFfldSfv8jsSa9bspQXzKSAGb2LUmlcj2XmbFtA79qZj8MPA94U77o81LgBjM7G7ghfx8EwRrHgK5VSm1rjJslfUDSefn2PjKBZyhDr9TM9prZV/PXB4BvkK1vupDseZ/8788cVdeDIDjumFJf718E7iDLvvhmsoXrbyxz4EhzlJLOAp4DfAk41cz2QmZMJZ2SOGYXmVsSTzntjFFOFwTBMWINGsGhmFkT+MN8G4nSY2dJG4GPAm8xs+JET7pzu81sh5nt2HLitlH7FwTBxBFm5bb1QqkRZT7h+VHgr8zsY3nxg5K256PJ7cBDw9rpWIV9zaKr4WMHi/Op1YQJ3zznT+rXE65sG+tFN6/FhHtk6m4sVvxsfK3ZYvlMxw+86gUQBug6gWQBV7hJ/aqp7QsXlprYHyXobsot0Tvfoh8sNxVYOFVfDUfcS/SvesR3L9WmomgIUHVENc8VdTnaCZFn4YRiea3qf75eAGGA+Xm/7WazmLHy8OHiZ1NP9G0UDOhO4YhyJQwdUUoS8AHgG2bWO2S9Bnht/vq1wCfG370gCCaOTa2Yc9SUudLnA68BXijplny7ALgcOF/SXWQ+k5evYj+DIJgg0yTmSKpKeoOk35L0/IF94wmzZmafJ70k76fKnCQIgrXE1M0//jkwT5aS9o8l/V8z+5V8389SIsza+hk7B0FQiiVf72kZUQLnmtmrzeyPyHzJN0r6WJ4zp9RFhKEMgqDAlKneT7owmVk7z+N9C/BPwMYyDUzU17vTEY8fKqrThxaKNzylene6vrqdqr845wTXTbhHptTwmnz3vlnHnXJzzf/ytOd9N7mZuh9wuL5YVM9rNV9RV8Jdr5JInVtNpab1SATd9Rv21WPN+O56SWrO55AKTlxPnDPl8lgvflG8tMMADSVcV+ub3PKTZ/cVu1fx/w+lVKBp/3ranWJ5c7H4v9BojGfsM8I3ZC1ws6SdZvappQIzu0zS/RyPeb2DIDj+MTRViraZ/ZdE+fuB95dpIwxlEAQF1tBjdWkk/axTvA+4zcyWXQcehjIIggJrSKgZhdcBPw78c/7+POCLwNMlXWZmf5k6MAxlEAT9WHpKeI3TBX7YzB4EkHQq2Rzlc4EbgTCUQRCUY8ypII4nzloykjkPAU83s8ckLZsSYqKGstmG7z5QLG+1ygeBTal6mzf6H+z+w44ymLjqrZv8tjc2/HvY6haVyI4lUt5WEsFe636/O04w2bmEFukpuQCNhLptM0W/YTr+NSZTyjpDjqSvd4qEMu1hKV/0hOKvZsKPfMbxsU7cp5R/+XzDuX9As1FU9+dr/v3bNOuv3mh3/PJO1/mebCiWjSlb7VTOUQKfk/RJ4G/z968AbpS0AXhiuQOnR9oKgmBMiI6V21Z0FmmrpOsl3ZX/PdGpc6akf5b0DUl3SHpzz753SLpvwLV6Od4E/AVwDlm4yKuAN5nZITN7wXIHhqEMgqAPY2ILzstkSUhlWFjiXWZ2Tr4tm1snz83zebKF5p8BbszLhhKGMgiCAmblthUyNEvCMhkWRkbSq8j8vV8BvAr4kqRXlDk2xJwgCAqMIOZsk3Rzz/vdZra75LGlsiQsMZBhYYlLJP1X4GaykedyCQ7fBvy7pTWTkk4mG1n+3bCOhqEMgqCf0ZYHPWJmO1I7JX2GLIPrIG8bpUuJDAvvBX6LbLbgt8jS0f63ZZqpDCwsf5SST9UTNZTtVpeHHy6qkQsLRfUz5Y580jbfR7bdTii/jqo8O+v/WlYSaWkXE23P1ouqaMr1a3PDv6BWIgJ2u+KkIp0pH20coNJNqMrOM5ObHheSaWk9X2ovMjng+24vh5OalcWE+t7w/dktEfVdzhdLtUT8gIS1mG34g5bNzkqFas2/f5V5//tQr25wy7vOfGC1UmyjXl25l7YBXU9lP5q2zF6U2iepVJaERIYFepf65BkVPzmkO5+SdB1ZLm+A/0zJnOExRxkEQYEuKrWtkKFZEpbJsEBuXJd4OXD7ciczs18DdgPPAp5NNk3w1jIdjUfvIAgKjEGoKcPlwEckvQ74LvBKAEmnAe83swv4foaF2yTdkh/367nC/XuSziEbBN8LvGHYCc3so2Sj05EIQxkEQR82oQjnZvYoTpYEM7sfuCB/ncywYGavKXMeSQfIjGlhV9aMbR7WRhjKIAj6mTJfbzPzA4iOwFBDKelM4ENkylWX7Ln+3ZLeAfwC8HBedWk4nKTV6vDAfUVPoYOPF1OONub8SfqFI/6k/myi/patRbey6gH/17K56LsfnrAxIRTViuXteb9uq5sob/jn3FgrBpNV4wS3brWbEHl83Yt5xzXPaokgvynXQU/MSdVN/Ncl3RI94UapyMwJkSch0JjX70TKW2v4/asdLgboBZhpFIWYSuIam6l0tQmXx5mqI3g6A61qIsj0qEzo0XvNUGZEubQy/quSNgFfkXR9vu9dZvbO1eteEASTxmDF7onTRpksjHuBpUWhByQd9cr4IAjWBjGi7Gek5UHOyvhLJN0q6UrPoT0/ZpekmyXd3Dz86Io6GwTBZJiQC+OaobShdFbGvxf4QbJIHHvJVsUXMLPdZrbDzHbMzJ80hi4HQbCamGUL3Mts64VSqre3Mv4oVsUHQbBGWE+jxTKUUb3dlfFLrkf526Gr4gFazRYP3LO3UH5o34FCWWM24dp3hu8332r67va4kpsAAA6uSURBVHrNM7cVyqqJ3LbNpi8TLyz4CurGDcV2Ok4wX4DZRMDhyqaEIuwomq2q/3HNasFvu+orrm0nRW6l46utlWpCPXZU8mpCPU6llE0F3XXdEhPqtqVcLBMquS0650zUVaLtSsrlsV38HFKphDfX/O/aTMN3JW1b8ZwzKq6MqKp8EOzl6ExZvtqVUmZE6a6MBy4edVV8EATHP0vxKIPvU0b1Tq2ML+VMHgTBGmOdCTVlCM+cIAgKTJNnzjgIQxkEQR/Zo/ex7sXxRRjKIAgKhKHsZ6KGstvpuAq3x+KCr/498O09bvnsJj/g6YNO2dzGVPBfXzFsLabqF5V5S6SrbSWCBTcSPtbNdrGdrQnf7Sa+Cqt6QlF3/NHrnUR610QE5Yqj5mrTFr+NRFraVEpZN+huQmlOqdvJvK1eX5r+qoGUrZCTlhagsa/4bbPEqoFUUOXFmaGBbJ6k6qxUqKaCNY+Cheo9SIwogyDoI4twfqx7cXwRhjIIggLx6N1PGMogCAqEoewnDGUQBH3YlAXuHQeRXCwIggJmVmpbCZK2Srpe0l3531QEsnsl3Sbplt4c4mWPHwdrbkSZ+nCO7D9YuryWUETbT93ulreavr/u4cNFhbHdnnfr7q/7v0kp3/AtG4sqeb3qq61K6LP1iq/idyuO6p2IuF2b8RXhRnN/oaxa8a/FEteuGf+cntLuRSaHhO82+Oo20D1U9EdX3Vemu03fv7xW9a/TvUov9S7paPC1ZsJf3sNLO5zynx+RhJv7uLkUuMHMLpd0af4+lRXxBWb2yAqOXxExogyCoI+ysSjHMI95IXBV/voq4GcmfHxpwlAGQVCga+U2YNtSYO582zXCaU5dikCW//VDg2Urlj4t6SsD7Zc9fsWsuUfvIAhWnxFGi4+Y2Y7UTkmfIUtMOMjbRujO883sfkmnANdL+qaZ3TjC8SsmDGUQBAVsTLK3mb0otU/Sg0txbSVtBx5KtHF//vchSR8HzgVuBEodPw7WpaFst/wJ7713f88tn9u80S2fd9wmD+73U8pu3Jxyg/RdLw8vFAWDI01f/Ngw63+pWwmh6ISZotvkBic9LkBl1k/fMVcttmEJMae66LdNwj1SjrtiKqVsMo1twi3RE24s4YaihBBjqRS5h4vuuUqIZKon+pe4h5VW8Zw68FixrJ3o2wjY5FwYrwFeC1ye//3EYAVJG4BKnthwA/Bi4LKyx4+LmKMMgqBAt2ulthVyOXC+pLuA8/P3SDpN0lK821OBz0v6GvBl4B/M7FPLHb8arMsRZRAEaSYVZs3MHgV+yim/H7ggf30P8OxRjl8NwlAGQdBPRDgvEIYyCIIBjG5Yyj7CUAZBUCChs61byqSrnSWT4mfy+n9nZm+XtBX4G+AssiyMrzKzx1evq8eOUdwjUwGHD5+4yS2vJJTVVqvortjp+K52h2b8Njpd/+OtqKjEtruJdLoz/n9Mp1H+N7Yhv27Kja/qiATWSKSlTfjapcZDnltiUt1OiBWdfUX3TYCqo553O37dyiE/gHVlNhGd2VG9O0eKqwncoMcjYgadTowoeymjejeBF5rZs4FzgJ2Snsf3/SzPBm7I3wdBMAVMIijGWmKoobSMpaFTPd+MCfpZBkEwOYyRXBjXBaXWUUqqSrqFbOX79Wb2JUr6WUrateQH2lrcN65+B0GwWlg27VBmWy+UMpRm1jGzc4AzgHMlPbPsCcxst5ntMLMd9YbvtRIEwfHFhKIHrRlGUr3N7AlJnwV2MkE/yyAIJssYvG6mijKq98lAKzeSc8CLgN9lgn6Wa4kDjz4xUvnBx331c9tp2wplD2/yA/du2eoHC953wFfJD55U9NPeNO/XPbLBL99YL6qw7YQSPlP3Ff/5ht/v2UZx8UTtsD9tU0mlsU2klPWC7qZ8t1PqdjLQr9dOQpVPrr5Z8H3AO4cOF8uOFOuOR/U2uqF691FmRLkduEpSlexR/SNm9klJ/w/4iKTXAd8FXrmK/QyCYILEgvN+hhpKM7sVeI5TPjE/yyAIJst6WvpThvDMCYKgD7OYoxwkDGUQBAViQNlPGMogCPowMzoTity7VghDeYwZRSU/8SlFJRzgUCKq+v4tvqr8xBNFf+ItW3wl96QTiwo5wGZHJV84wf86nTzrK9bNhDK92YmePtPwI8FX24l0uvsedMvdhcNOZHLwfbchoW6Dq3B3E9H0SZRbQiVvOXEFWgeKSni3PZ48s+tpMXkZwlAGQVAgDGU/YSiDIOhnnflxlyFy5gRB0IcxGV9vSVslXS/prvzviU6dH5J0S8+2X9Jb8n3vkHRfz74LVtShZQhDGQTBAOVCrI1hreXQUI1mdqeZnZPHmvgx4DDw8Z4q71rab2bXDh4/LuLRew3x+AOPuOULh31BY/9jvpjTPG1roazV8sUS8AWXhWbRFbBW9VOz1iuJdL81XxSp1oqCRCUR5FedhChS9cUpnCC9qZSyqaC7SbdER6BZeND/zFTxxygpMaa5r5iu98DeouDXXmi5x4+EMSnV+0LgvPz1VcBngbcuU/+ngG+b2b+tbreKxIgyCII+JvXoTclQjT1cBHx4oOwSSbdKutJ7dB8XYSiDIOhntHiU25bizebbrt6mJH1G0u3OduEoXZLUAH4a+Nue4vcCP0iWeWEv8Acruu5liEfvIAgGGCkL4yNmtiPZktmLUvskjRKq8SXAV83syQWyva8lvQ/4ZNlOj0qMKIMgKDChR++lUI0wPFTjxQw8dufGdYmXA7evtEMpwlAGQdCHMbHkYpcD50u6Czg/f4+k0yQ9qWBLms/3f2zg+N+TdJukW4EXAL+80g6liEfvKWCUdLoAzUPFNKf7H/PdIJsLJ7nlJ2wpquH1hHqcxXsusmnWV6Yr80XFtZloe3PNb7vS9dVfL0Wu6v6qgVRK2aQe7KjeKXXbEu6Rhx/23T337SkGMz7yePFz7LbG4MJo0GmvvuqdCtVoZvcDF/S8PwwUvohm9ppV7WAPYSiDIBhgfaWiLUMYyiAI+jBLj3jXK2EogyAoEIF7+wlDGQRBgXj07meo6i1pVtKXJX1N0h2SfjMvn5hDehAEk8PM6La7pbb1QpkRZRN4oZkdlFQHPi/pH/N97zKzd65e94LV4NC+opp75GAxCCzA3Abf19tbQzc/73+dvBSxAO2Or3rXq0W/85Rf+Eyj6ZYvzmz2+9Is+kyr4vevMusr6qmUsl7Q3ZTvdkrdfuK7j/nl3ynWP3Rv8drbCT/0Uena+jGCZSiThdGApXUm9XyLcXkQTCsWgXsHKbXgXFJV0i1kLkbXm9mX8l1DHdIl7VryA20t+r+iQRAcPxjlvHLWkzEtZSjNrJPHgzsDOFfSMynpkG5mu81sh5ntqDf8Rc1BEBxfTMgzZ80wkgujmT1BFjNup5k9mBvQLvA+4NxV6F8QBJPGoNvtltrWC2VU75MlbclfzwEvAr45SYf0IAgmh2F0O51S23qhjOq9HbhKUpXMsH7EzD4p6S8lnUMm7NwLvGH1uhmsNqkv/ZFDvsIrJ1J4s+lHVG93fFW50y22AdC1YvlM1ffdblsikvkIVFqJ9LOJ8s4hf4WAl1LWi0wOvu82+Oo2+Ar3qhFiToEyqvetwHOc8ok5pAdBMFnCUPYTnjlBEAxgsY5ygDCUQRD0YfHoXSAMZRAEBSJ6UD9hKINlaSZS4TZmG4Wyw4f91LHNxYTgssEXc6qV4j+p4dedUTF4LUC1kxBonLV/OuC7DXaO+G13jvj3pHWgKPJ4KWXBD7oLExZtUpitK0W7DGEogyDow4gwa4OEoQyCoJ8I3FsgkosFQTDAZHy9Jb0yD93YlZRMeStpp6Q7Jd0t6dKe8q2Srpd0V/7XjTcxDsJQBkFQwKxbalshtwM/C9yYqpA7uryHLK/3M4CLJT0j330pcIOZnQ3ckL9fFeLROwiCPrLAvasv5pjZNwAkX6jLORe428zuyeteDVwIfD3/e15e7yqyOBRvXY2+TtRQHtr3rUe+8Mnz/i1/uw14ZJLnPwbENU4Ha+kan7rSBg7t+9Z1//LJ87aVrD4r6eae97vNbPdK+9DD6cD3et7vAZ6bvz7VzPYCmNleSaeM8bx9TNRQmtnJS68l3WxmyXmJaSCucTpYD9fYi5ntHFdbkj4DPMXZ9TYz+0SZJpyyiUvy8egdBMGqYWYvWmETe4Aze96fAdyfv35Q0vZ8NLmdLLD4qhBiThAExzM3AWdLepqkBnARcE2+7xrgtfnr1wJlRqhHxbE0lOOcxzheiWucDtbDNU4cSS+XtAf4ceAfJF2Xl58m6VoAM2sDlwDXAd8gC/N4R97E5cD5ku4Czs/fr05f11M49yAIgqMhHr2DIAiGEIYyCIJgCBM3lCl3pLVOnrL3IUm395RNzMVqEkg6U9I/S/pG7nr25rx8aq5T0qykL0v6Wn6Nv5mXT801BqMzUUM5xB1prfNBYHD92cRcrCZEG/hVM/th4HnAm/LPb5quswm80MyeTZaKeaek5zFd1xiMyKRHlE+6I5nZIrDkjrTmMbMbgcHAhheSuVaR//2ZiXZqzJjZXjP7av76AJkKeTpTdJ2WsZQlrJ5vxhRdYzA6kzaUnjvS6RPuwyTpc7ECVs3FatJIOoss6dyXmLLrlFSVdAvZAubrzWzqrjEYjUkbyuPCHSlYGZI2Ah8F3mJm+491f8aNmXXM7BwyL5BzJT3zWPcpOLZM2lAu5440jTyYu1ax2i5Wk0JSncxI/pWZfSwvnrrrBDCzJ8gi0uxkSq8xKMekDeVy7kjTyMRcrCaBsnhYHwC+YWZ/2LNraq5T0smStuSv54AXAd9kiq4xGJ2Je+ZIugD4I6AKXGlmvz3RDqwSkj5MFhtvG/Ag8Hbg74GPAD8AfBd4pZn5mazWAJJ+AvgccBuwFLX118nmKafiOiU9i0ysqZINJD5iZpdJOokpucZgdMKFMQiCYAjhmRMEQTCEMJRBEARDCEMZBEEwhDCUQRAEQwhDGQRBMIQwlEEQBEMIQxkEQTCE/w/zB26gl7OH6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(\n",
    "    np.log2(S / C)[1:-1, 1:-1],\n",
    "    cmap='coolwarm',\n",
    "    vmin=-1,\n",
    "    vmax=1,\n",
    "    \n",
    ")\n",
    "plt.colorbar(label='log2 obs/exp') "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "_draft": {
   "nbviewer_url": "https://gist.github.com/9abdc3597f132d9b24889ab83bf70fbb"
  },
  "gist": {
   "data": {
    "description": "saddle_plot_tutorial.ipynb",
    "public": false
   },
   "id": "9abdc3597f132d9b24889ab83bf70fbb"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  },
  "nav_menu": {},
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
